FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 17, N o 2, 2019, pp. 217 - 242 https://doi.org/10.22190/FUME190227027F © 2019 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper * METHODS OF THE PATTERN FORMATION IN NUMERICAL MODELING OF BIOLOGICAL PROBLEMS Alexander E. Filippov 1,2 , Stanislav N. Gorb 1 1 Functional Morphology and Biomechanics, Zoological Institute, Kiel University, Germany 2 Donetsk Institute for Physics and Engineering, National Academy of Sciences of Ukraine Abstract. Evolution of different systems can be described in terms of their relaxation to the minimums of some effective potential relief. This observation leads us to face us with a question how to generate corresponding potential patterns which describe adequately various physical and biological systems. In this review, we present a number of different ways of generating such potentials demanded by the problems of different kinds. For example, we reproduce such a generation in the framework of a simple theory of phase transitions, automatic blocking of the growing phase nucleation and universal large scale structure. Being frozen at late stages of their evolution they form majority of meta-stable structures which we observe in real world. Counting on above-mentioned universality of naturally-generated fractal structures and their further utilization in numerical simulations of biological problems, we reproduce also formal algorithms of generation of such structures based on random deposition technique and Fourier-transform approaches. Key Words: Pattern formation, Phase transitions, Large river effect, Nucleation, Biological applications, Frozen kinetics 1. INTRODUCTION A large problem of experimental biology is that we are dealing with systems which are very fragile and sensitive and behave statistically, which means that often no kind of experiment is possible, and if so, only at extremely high costs and usually without proper control samples/experiments. This is one of the most important fundamental differences between experimental biology and other experimental disciplines dealing with non- biological matter. Again, the most helpful approach in this context would be modeling, which can then be used as an experimental platform “in silico” in order to predict Received February 27, 2019 / Accepted May 18, 2019 Corresponding author: Alexander E. Filippov Donetsk Institute for Physics and Engineering, National Academy of Sciences of Ukraine, Donetsk, Ukraine E-mail: filippov_ae@yahoo.com 218 A.E. FILIPPOV, S.N. GORB outcomes in real biological experiments. In some cases, a „virtual‟ experiment may even be the only viable possibility. These considerations lead to the insight that in many instances experimental biology can profit from new methods made possible by novel modeling techniques that take advantage of the recent developments in modeling and visualization. In many concrete examples we will need numerical generation of the effective potential. If effective potential is known and dynamic equations are properly written, the system is attracted to the correct configurations by itself. The potential can appear in real space, or in some imaginary space of parameters. In some sense, it does not matter. However, it can matter and it should matter from either physical or biological points of view for particular problems, where the potential has some specific meaning. In this review, we will present some practically useful and technically convenient examples of effective potentials (or patterns, if we are talking more definitely about their formation in 2D or 3D spaces). From the very beginning, one can divide them into two large classes. 2. SIMPLE THEORY OF PHASE TRANSITIONS AND PATTERN FORMATION One kind the potentials are abstract numerically generated formations which are simple enough and convenient for the theoretical numerical simulations. The other ones are provoked by the studies of the specific, experimentally observed systems. In this latter case, the degree of simplicity of the potential or the procedure of its generation is not so important because in such a case we are more interested in the potential (surface, pattern, density distribution, etc.) by itself at the final stage of the study, or maybe even as a goal of the particular research. For example, if we would like to know how this particular picture could appear self-consistently. Keeping this in mind, we will continuously jump between these two opposite limits. Sometimes, using very abstract models to illustrate general ideas, or simply to generate the equations, reproducing what we have observed, even without complete motivation. But, sometimes, we will devote almost the whole study to an extraction of the equations or procedures from particular nature of the problem under consideration. The same note can be done about dimensionality of the surface (potential, line, etc.) generated by a procedure. In some cases it is enough to create 1D line with specific properties. It is important to note that even for such “trivial” case our motivation can be different. For example, it can be made just for simplicity of further application of the generated object. However, in many cases, it can be even more correct than to generate and apply 2D pattern and use it. It can be proven that for many problems of contact mechanics, real 3D problem can be exactly transformed to the 1D model! It is so-called Method of Dimensionality Reduction (MDR) invented and actively developed last during decades by Valentin Popov and co- workers [1]. The cases of the MDR application we will mention in this paper and specially discuss them in corresponding sections. It should be additionally mentioned that the absolute majority of patterns, we have used in studies of biological problems under consideration, were generated by simulations of physical or chemical kinetic processes. Below, we basically report on the methods of physical kinetics. For example, the ideas from fluctuation theory, phase separation, and phase transitions can be applied for pattern generation. Of course, it is practically impossible Methods of the Pattern Formation in Numerical Modeling of Biological Problems 219 to reproduce all the mathematical foundations of this branch of science, which are published in hundreds of the handbooks starting from the 19 th century. Thus, in many cases, we will simply suppose that it is either known or can be proven. Independently of the accuracy and correctness of the terminology, one can start from a very formal and methodically transparent approach. Let us suppose that there is fluctuating density ρ(x) randomly distributed originally at t=0 in one-dimensional space {x}. Random initial distribution means that fluctuations δρ(x,t) on average are equal to zero <δρ(x,t)>=0 and originally independent at every point and time of the system <δρ(x,t)δρ(x´,t´)>=Dδ(x-x´)ρ(tt´). It is so-called δ-correlated noise. In the majority of real systems the densities of fluctuations interact. These interactions can cause mutual attraction and enforcement of the fluctuations leading to the nucleation of the regions with nonzero density ρ(x) and their expansion in space. One of the most accepted theories of this process is Landau theory [2] of phase transitions (which can be treated as a combination of main interactions within the system in the framework of the simplified “mean field” theory). Main idea of Landau theory is that energy of the system can be expressed as a function of the so-called “order parameter” density. In particular, it can be a functional U(ρ(x)) of spatially distributed real density ρ(x). The simplest variant of the theory corresponds to the case, when local form U(ρ(x))→U(ρ) of the energy (depending on uniform variable ρ) has two minimums corresponding to “ordered” ρ=ρ0≠0 and “disordered” ρ=0 states. Depending of the relationship between the energy minimums U(ρ0) and U(ρ=0) one or another state is energetically preferable. For example, if U(ρ0)=Dδ(xx´)ρ(tt´) tend to the ordered state with <δρ(x,t)>=0 and <ρ(x,t→∞)>=ρ0=const. In this place, the reader can ask us: why do we need so complex description, if after all everything tends to the simple constant ρ0=const? The main problem here lies in the words “after all”. We are not interested in trivial final of the kinetic process, but in its intermediate state. Absolute majority of real surfaces and substances are “frozen” in an intermediate states of the kinetic process and namely at these stages they produce Methods of the Pattern Formation in Numerical Modeling of Biological Problems 221 practically important patterns with non-uniform ρ(x) with nonzero gradients ∇ρ. Moreover, in many cases (for more complex systems than we study right now) even final stage of the kinetic process are domains and other structures with non-uniform distributions ρ(x). Below we will study some of them, but now let us return to the simplest (!) case. Combining Eqs. (1) and (3) one can write energy functional in the following form [3]: 2 2 3 4 2 [ ( )] [ ( ( )) ( ) ( ) ( ) ...] 2 2 3 4 [ ( ( )) ( ( ))] 2 c b g x dx x x x x c dx x F x                      (4) However, if it is still not enough to write down complete equation for the fluctuating density [4]: ( , ) / [ ( , )] / ( , )x t t x t x t        (5) If ρ(x,t) really fluctuates, it is not only random value <δρ(x,t)δρ(x´,t´)>=Dδ(x-x´)ρ(t- t´) with intensity of D at the beginning of the process, but it continues to fluctuate further. It means that Eq. (5) should include also a source of random fluctuations δ(x,t). Mathematically, it has the same properties as the initial condition for the fluctuations of density: ( , ) 0x t  , ( , ) ( ', ') ( ') ( ')x t x t D x x t t       . (6) This random source stays in the right hand side of the equation of motion ( , ) / [ ( , )] / ( , ) ( , )x t t x t x t x t          (7) and cause generation of new fluctuations during the whole process of the density nucleation. Eq. (7) is already equation of motion for this simplest (!) model. Taking into account explicit form Eq. (4) of the functional Φ[ρ(x)] and performing the variation δΦ[ρ(x,t)]/δρ(x,t) one can write the form of the kinetic equation that is more convenient for the solution: 2 3 4 ( , ) / [ ( , ) ( , ) ( , ) ( , ) ] ( , )x t t c x t x t b x t g x t x t               (8) where Δρ(x) is Laplasian operator which in 1D case is equal simply to second derivative Δρ(x)=∂ 2 ρ(x)/∂x 2 . Starting from the random distribution, the density evolves with time according to the Eq. (8). The barrier in the local part of energy F(ρ(x)) does not allow majority of the density fluctuations to grow immediately to the minimum corresponding to the nonzero equilibrium density ρ0=const. The only relatively large fluctuations can pass the barrier and grow. They grow in both senses: their amplitude at maximum tends to ρ0 and they expand in space [5, 6]. The same logics can be easily expanded into 2D space where corresponding pattern of the density will appear. Typical intermediate configuration with different nuclei of nonzero density is shown in Fig.2. If there is physical reason to freeze this configuration, we would get almost static representation of the non-uniform pattern of density. The simplest possible reason for this may be very small damping constant  which controls 222 A.E. FILIPPOV, S.N. GORB rate of the kinetic process in Eq. (8). If the characteristic time 1/γ of evolution of the distribution ρ(x,t) to the equilibrium ρ0 is extremely long in comparison with other times of a problem, one can ignore the evolution of the density at all and treat its intermediate distribution as practically static. It has some ideological analogy with the galaxies developing for billions of years. It is much longer than our life and majority of the processes on our planet. So, the galactic patterns can be treated as practically static ones. Fig. 2 Typical intermediate pattern appearing as a result of solution of Eq. (8) However, in many cases the kinetic process may be really stopped by some natural reasons at an intermediate stage. Normally it happens when the system is actually open to the external influences or consists of a number of interacting subsystems. In the first case, external forces can cause change of the coefficients of energy expansion Eq. (4) already in the course of the density evolution. As result, the system starts its development to some equilibrium at the given form of F(ρ) and finishes in absolutely another equilibrium corresponding to modified functional F(ρ). In the second case, alternating interactions between the subsystems can lead to the growth of different densities in different domains of the space. When such growing domains meet each other, they mutually block further expansion and build a static (or maybe very slowly shifting) domain wall between them. Such blocking leads to the fixation of the domain structure and on large scale forms a pattern including many domains of different sizes. The simplest way of obtaining such domain structure is to use expansion F(ρ) with even terms only: 2 2 4 6 2 [ ( )] [ ( ( )) ( ) ( ) ( ) ...] 2 2 4 6 [ ( ( )) ( ( ))] 2 c g u x dx x x x x c dx x F x                      (9) Corresponding equation of motion with the functional Eq. (9) leads to the typical domain structure shown in Fig. 2. Methods of the Pattern Formation in Numerical Modeling of Biological Problems 223 Fig. 3 Formation of the domain structure in the kinetics. Two different kinds of domains are shown by red and blue colors, respectively 3. AUTOMATIC BLOCKING OF THE NUCLEATION AND FREEZING OF THE PROCESS In both previous cases, telling about almost frozen kinetics, we have supposed that it this process is simply much slower than other ones in a particular model. However, it was shown more than 20 years ago that since the formation of new phase nuclei includes processes that prevent their appearance and growth in other regions in space, it should result in autostabilization of an intermediate mixed state. One can call it “automatic blocking of the nucleation” [7]. Normally it is caused by the nonlocal (long-range) interactions which are caused in an ordered process by itself. From the introduction, we remember that numerical modeling does not like nonlocal interactions because of an extreme time consumption required for these calculations. However, it is possible to show that the effect of blocking long-range interaction can be approximately included by some local additions to the above models. Various mechanisms for the formation of an effective long-range interaction in such systems can be analyzed. Here we will mention only a couple of them. Normally, formation of the substance (and the surface, which is needed here) includes the reaction of a crystal lattice (striction) to a change in magnitude of order parameter ρ(x,t) during a phase transition. One can show that in a relatively simple case of an isotropic medium and quadratic striction, the local energy functional of the free energy is modified in the following manner [3]: 2 2 2 2 2 1 1 , [ ( )] [ ( ( )) ( ( )) 2 1 ( ( ) ( ) )] 2 3 u ii ii k k ll l k c x dx x F x k g x u u u u                 (10) If the lattice vibrations manage to follow the variations of ρ(x,t), we can utilize condition δΦ[ρ,u]/δuik=0 to eliminate variables uik and get, after some standard mathematical transformations, an effective functional solely in terms of field ρ(x,t): 2 2 2 [ ( )] [ ( ( )) ( ( )) ( ) ' ( ') ] 2 2 c x dx c x F x x dx x V            (11) Here ( ( ))F x is the renormalized local form of F(ρ(x)) with the same structure as the original function F(p) (we shall henceforth omit the tilde), and constant κ is defined by 224 A.E. FILIPPOV, S.N. GORB expression 2 1 1 ( / 2 )[( / 2 2 / 3) ( / 2 2 / 3) ]q V k k P        . Constant P is determined by the external pressure or other constraints which prevent free expansion of the crystal (twins, defects, etc.) and, in turn, fix the sign of κ. When P>μ, κ>0; otherwise κ<0. Nonlocal construction ∫dx[ρ(x) 2 ∫dx´[ρ(x´) 2 ] in functional Φ[ρ(x)] generates a term with a long-range effect in the equation of motion for the field variable 2 ( , ) / [ / ( , ) ' ( ', ) ] ( , ) 2 x t t c F x t dx x t x t V                   (12) whose presence significantly accelerates or slows down (or even totally stops) the ordering process, depending on both the magnitude and sign of κ. Typical picture of blocked almost static domains is shown in Fig. 4. Fig. 4 Formation of the domain structure in the kinetics at the presence of self-blocking. Two different kinds of domains are shown by red and blue colors, respectively. It is important to note that the last stage here is the practically final one and does not develop further with time. Let us discuss on one more example of an interaction, which leads to a similar model. The local variation of the order parameter is accompanied by the evolution or absorption of heat (depending on whether the transition is to the low- or high-temperature phase). This results from heat conduction in heating (cooling) of the surrounding regions of space, which, of course, slows down the transition process in all cases. This mechanism seems to be universal, and its effectiveness is determined only by the relationship between the rates of the nucleation and heat-conduction processes. The local heating (cooling) of a system in a region, where a nucleus appears, can be taken into account by assuming that quantity τ in expressions (9) is a function of position and time. The kinetic equation for the order parameter should be supplemented by an equation which describes the evolution of τ=τ(x,t). The latter equation should be a heat-conduction equation with the heat removal and with a source β[ρ], whose intensity is proportional to the rate of change of the free energy, i.e. β[ρ]~∂ρ/∂t·δΦ/δρ. As a result, we have the following system of the connected equations [7]: ( , ) / [ ( , )] / ( , ) ( , )x t t x t x t x t          ; (13) ( , ) / ( , ) ( , ) / [ ( , )] / ( , ) ( , )x t t x t x t t x t x t x t                (14) Methods of the Pattern Formation in Numerical Modeling of Biological Problems 225 Before starting discussion of the results of the numerical experiments, we show that with some roughening of the model, the mechanism under consideration can be described in terms of a single field variable ρ(x,t), which evolves in accordance with an equation similar to Eq. (12). The physical arguments, which lead to a functional like Eq. (11) in this case, too, are fairly simple. Each growing domain of the new phase creates a non- uniform temperature field τ(x,t) around itself. Owing to the heat conduction, the temperature at other points in space deviates from trial temperature altering the conditions for the growth of other domains at those points. This signifies the appearance of an effective long-range field accompanying the nucleation process in the system. Relating the variation of the temperature field to the order-parameter field ρ(x,t), we obtain an energy functional Φ[ρ] like that in Eq. (11). In fact, when the fluctuations of τ(x,t) are "turned on" in a system with a temperature equal to the heat-bath temperature τ0, after a unit of time the mean value of τ(x,t) deviates from τ0 by 1 2 0 0 1 1 /dx t V V            (15) The mutual influence of the domains of the new phase becomes significant, when they become so large (and this is seen from the results of numerical experiments) that the energy of the domain boundaries between the ordered and unordered phases can be neglected. As result, we arrive at a functional like Eq. (1), which was obtained to describe striction effects in the kinetics of a first-order phase transition. In some cases, this makes it possible, in principle, to disregard the specific mechanism for realizing the long-range effect accompanying the first-order phase transition and to formally analyze models with nonlocalities of the general form. 4. LARGE SCALE STRUCTURE OF THE FLUCTUATING FIELD. UNIVERSALITY AND SCALING Trial structure of the energy functional has direct relation to the microscopic interactions in the system. In many cases, it can be even analytically derived from the microscopic theory. It means that the coefficients of the expansions used in the kinetic equations above have well defined specific values and in turn should completely determine density distributions ρ(x,t) and as result the structure of the contact surface. However, it is well known that absolute majority of real surfaces have practically universal (scaling) structure with the power low distribution of the relief. It means that if there is no special reason to produce another structure, it should appear due to universal kinetic process which makes difference in the initial energy functional negligible. As we saw in the previous section, ordering transition is anticipated by nonlinear excitations, which can be interpreted as nucleation centers. The kinetics of the first-order phase transition in different physical systems has been the subject of intensive studies. As a rule, the ordering of a metastable disordered phase is due to the fluctuation production and finally to the growth of the nucleus of the stable phase. In a first-order phase transition, there is a change in some order parameter between these two phases, which lowers the free energy as the new phase forms. 226 A.E. FILIPPOV, S.N. GORB The corresponding local energy density F(ρ) must have a metastable minimum and be energetically favorable. However, the free energy is transformed due to the fluctuations. This change is especially essential in the critical region at a second-order transition. It is well known from the theory of critical phenomena that the fluctuations manifest themselves by renormalization of critical exponents. But, the renormalization group (RG) method allows one not only to perform purely numerical calculations of critical exponents, but also to predict some qualitatively new effects which could not be obtained within conventional approaches, e.g., within the Landau approximation, applied above. Among them, there are qualitative effects, such as the fluctuation-induced first-order phase transition. This effect takes place in some anisotropic systems, where the renormalized free energy F(ρ) undergoes transformations, which are typical for the first-order phase transition. The same kind of the nucleation centers can be found in a fluctuation-induced first order transition. However, one can expect that even in the situation, when the fluctuations are not strong enough to change the transition order, they manifest themselves somehow. The mean field model is very convenient for analytical study, but it reduces the fluctuation interaction, and it is impossible to control a correction which should be done to the free energy with account of the neglected fluctuations [8]. A more correct approach was developed almost 50 years ago. It states that the fluctuations renormalize energy functional Φ[ρ(x,t)] according to so-called renormalization group (RG) equation ˆ/ l R    . Here we do not need the complex exact form of the RG equation and use its simplified approximate version: 2 ( )ˆ [( 2) ( )] ' 2 ( ) ( ) ( ') ( ) ( ') d d r r R d r d r r d r r r r r r                                   (16) where d is the space dimensionality. The first term in Eq. (16) corresponds to a simple scale transformation of the density distribution ρ(x,t) and the second term appears due to the integration over internal fluctuations inside small regions after scale transformations. According to the general RG hypothesis in critical point of the phase transition (exactly, where the ordering of ρ(x,t) takes place) the functional Φ→Φ* tends to the fixed point * *ˆ/ 0l R     . Let us study now the time-space evolution of fluctuating density ρ(x,t) in this state. We plan to show that ρ(x,t) produces a well-pronounced large-scale structure in spite of its scale invariance on average. As always kinetic equation can be written in the form: ∂ρ/∂t=−γδΦ*[ρ]/δρ+δ. At every time moment, average probability W= to find density distribution ρ(x,t) is determined by functional w[ρ]=exp(−Φ[ρ]). This probability develops with time according to the equation: [ ] / / d w W t d r t          . (17) By combining Eq. (16) with condition *ˆ 0R  and ∂ρ/∂t=−γδΦ*[ρ]/δρ+δ, one can obtain that in the critical point time the evolution of probability W= is reduced to the simple scale transformation: * / [ ] [( 2) ] 2 d r W t w d r d r              (18) Methods of the Pattern Formation in Numerical Modeling of Biological Problems 227 Physically it means the following: at the point of the transition energy, the functional becomes the universal one Φ→Φ*. The kinetic equation based on this functional produces with time new realizations of ρ(r,t) with the same statistical properties, but larger scales. In the limit t→∞, the structure becomes scale invariant. This structure is generally similar to that found in the intermediate stage of the nucleation process at first- order phase transitions, but it never finishes in static ordering. To calculate practically some particular realization of the density in the critical point, one can solve numerically equation ∂ρ/∂t=−γδΦ*[ρ]/δρ+δ with the energy density found from the RG equation in its local approximation: 2 [ ( )] [ ( ) ( )] 2 c x dx F      , (19) where ( )F  is normalized to the critical temperature * 2 ( ) ( 0)F F F       solution of the local version of RG equation *ˆ 0R  : 22 2 2ˆ 0 2 d F F F RF dF                   (20) Despite of its visual simplicity, Eq. (20) is nonlinear and cannot be solved analytically. But, in contrast to functional equation *ˆ 0R  it is already an “ordinary” differential equation and it is easy to find physical branch F * of its solution *ˆ 0RF  numerically with a very high accuracy. Discrete data array F * (ρk) has big enough number N≫1 of points k=1,2…N, where each value of {ρk} defines unique value of {f * k}. With a good accuracy, it can be used in the kinetic equation: * / [ / ]t c f            (21) instead of the analytic formulae for the energy, otherwise normally used before. Repeating the simulation with random initial conditions and with the random time- dependent noise <δ(r,t)>=0, <δ(r,t)δ(r´,t´)>=Dζδ(r-r´)δ(t-t´) after long-time runs, we can get unlimited number of the realizations ρ(r,t). Direct observations of the simulation results show that every instant density distribution ρ(r,t) contains many nuclei of different size. Then, the wider particular maximums of density are formed, the longer time they survive in the general landscape. It means that at large scales the total process of transformation becomes slower and slower. Besides, stronger correlation between the densities in different spatial positions appears. In an extremely long run limit (when t→∞), density distribution ρ(r,t) tends to the expected scale invariant structure. One can calculate correlation function G(rr´) =< ρ(r,t)ρ(r´,t´)> at the fixed time moment t and find this scaling correlation function. It means that G(rr´)=<ρ(r,t)ρ(r´,t´)~1/|rr´| β depends as a power function on distance |rr´| between the points. This quite general result favors to the common observation that in many cases, irrespective of the specific features of the system under investigation, the evolution of the surface layer proceeds in a fairly universal manner. 228 A.E. FILIPPOV, S.N. GORB 5. CHEMICAL APPEARANCE OF FRACTAL SURFACES One more example of automatic generation of fractal surfaces comes from the case, where in the immediate vicinity of the smooth flat interface of the two media in contact, a dense layer of one or more reaction products emerges, as a result of a chemical reaction. It was observed that in the process of growing this layer becomes more and more porous and rough. Gradually, an essentially inhomogeneous but, as a rule, scale invariant structure is formed, and the laws governing the growth of this structure are characterized by fractal dimension and growth exponent [9]. In addition to arousing purely scientific interest, the study of corrosion-front growth attracts attention because of its importance from the standpoint of practice, since in some applications the problem is closely linked to that of raising the efficiency of electric batteries. For instance, when a lithium anode is placed into an electrolyte containing SOCl2 as an additive, due to the exceptionally high reactivity of lithium, a porous two- component layer of LiCl and SO2 is formed at the surface of the anode. The presence of such a layer leads to what is known as a lag effect, when the element is stored for a long time. Micrographs of the surface layer show that the layer can be considered as a combination of a relatively dense initial layer with a subsequent transition to a fractal structure with an ever increasing porosity. The highly universal properties manifested by different systems suggest that one can use universal growth models based on a combination of the ideas of continuum field theory and kinetic equations with a random source. Being fairly common in the theory of phase separation and fluctuation phenomena in phase transitions, the kinetic equations with a source of noise should be used cautiously in describing front growth, the reason being that, in contrast to phase transitions, where generation of the order parameter occurs in the bulk of the system, a random source cannot be considered additive. The generation of a finite density of components forming the front occurs only in the immediate vicinity of the already existing boundary. This means that in the case at hand the corresponding source in the equation must be multiplicative i.e., at least contain the density as a factor. However, in recent publications, devoted to theoretical studies of phase diagrams and transitions in systems with multiplicative noise, it was noted that the presence of such strong noise can have a dramatic effect on the ordered structure and on the phase diagram, and may lead to the emergence of new nontrivial phases. In our case, this means that the model equation should be written in such a way so as to exclude additional difficulties associated with this noise. From an experimental standpoint, the study of fractal corrosion structures is convenient since the corrosion front is observed directly in micrographs and the corresponding two dimensional distributions of density can be studied explicitly. At the same time, the processes involved are very complex, and notwithstanding the continuing efforts, the theoretical models still remain extremely simple, although they presuppose a numerical analysis of the kinetic equations. Usually, only the density of a single distributed quantity that is considered the most important in each specific case is involved. This is generally not the case in physicochemical processes, since usually two or more components participate in the reactions. No matter how subtle the description of a system by the single-component approach is, it is sufficient to replace the study of the system by an analysis of purely theoretical models. Given contemporary computer modeling Methods of the Pattern Formation in Numerical Modeling of Biological Problems 229 techniques, any attempt to reduce the problem to a single equation is more a tribute to the analytic tradition than a real necessity. Here we will demonstrate the feasibility of moving in this direction by the example of a two-component model formulated for the description of growth and corrosion of a broad class of porous surface layers initiated by chemical reactions. Below we examine the simplest two-component case, assuming, as an example, that we are dealing with chemical reactions that proceed in a system with a contaminated lithium anode. The complete picture of the reactions in such a system is fairly complicated and can be expressed as follows: Li→Li + +e - , 4Li + +4e - +2SOCl2→4LiCl+SO2+S. Actually, we are interested only in the formation of a front consisting of lithium chloride LiCl contaminated by the reaction product SO2 that concentrates near the surface. Bearing all this in mind, we can interpret the bare equation 9 as an initial equation for the evolution of the density of LiCl which we denote by ρ1(r,t). The corresponding coefficients and the source of noise will be labeled by the index “1”. By ρ2(r,t) we denote the density of SO2. We model the local repulsion of the reaction of the products LiCl and SO2 by a fixed-sign additional term in the effective energy of the system, V12(ρ1,ρ2), which in the lowest order that can be written V12(ρ1,ρ2)=Bρ 2 1ρ 2 2/2. Equation of motion for the first density becomes ∂ρ/∂t=−γ[c1Δρ1+ρ1(1−ρ1)][C1+δ1(r,t)]- Bρ1ρ 2 2. This equation must be augmented with an equation describing the evolution of the second component, ρ2. The second component ρ2, just like the first one, is generated as a result of the same reactions near the free (not contaminated by SO2) LiCl surface. This means that for ρ2, we must use the same generating term ρ1(1−ρ1) as for ρ1: 2 1 2 2 1 1 2 2 2 1 2 / [ (1 )][ ( , )] ( ).t c C r t B f                  Here we have allowed for the fact that although both densities, ρ1 and ρ2, emerge as a result of the same reaction, the rate of formation of the dense components in r ρ1 and ρ2 may differ, so that generally C1≠C2. Obviously, the terms linear in ρ2 cannot ensure that the increase in ρ2 is stopped and is stabilized ρ2→1 in the static limit. We must also bear in mind that far from the front, there is no spontaneous generation of ρ2, and hence the effective energy V2(ρ 2 2)=Bρ 2 2(1−ρ 2 2), whose variation yields the combination proportional to the function f(ρ2)=ρ2(0.5−ρ2) (1−ρ2). It contains a barrier that separates the two similar minima at ρ2=0 and ρ2=1. It can be shown that in the continuum approximation, with only the lowest harmonics in the energy (and hence Laplasian terms Δρ1,2 written in the equations), one has to add generation terms with the step-function cut-off factor | '| ( ' ( ') ) r r dr r a         . It turns on the generation, when the density in some neighborhood |r-r´|<σ exceeds critical threshold a. In other words, Θ→1 when ' ' ( ') r r dr r a      and Θ→0 in the opposite limit. Finally kinetic equations take the form: 2 1 1 1 1 1 1 1 1 2 / [ (1 )][ ( , )]t c C r t B               (22) 2 1 2 2 1 1 2 2 2 1 2 2 2 / [ (1 )][ ( , )] (0.5 )(1 )t c C r t B                     . In purely local approximation, the system Eq. (22) can be reduced to the form of ordinary differential equations: 230 A.E. FILIPPOV, S.N. GORB 2 1 1 1 1 1 1 2 / (1 )[ ( , )]t C r t B           (23) 2 1 1 1 2 2 2 1 2 2 2 / (1 )[ ( , )] (0.5 )(1 )t C r t B                 . Here the scale of time is normalized to γ=1. The global structure of the phase portrait is depicted in Fig. 3 which presents a physically interesting realization of the portrait. The local system of equations Eq. (23) actually describes the evolution of the densities at each point in space without allowing for interaction between different points. In this approximation, the interaction is taken into account only via the initial conditions. Specifically, as the front arrives at a point in space, both densities ρ1 and ρ2 begin to be generated at that point, so that the physical scenario in the phase portrait in Fig. 5 corresponds to the trajectories that emerge in the neighborhood of the trivial point ρ1=ρ2=0. Fig. 5 Phase portrait of local version of the kinetic equations describing growth of porous surface. Fixed points corresponding to possible combinations of local densities are marked by red color. Two well pronounced “large rivers” leading to the static final configurations are shown by the blue lines. Cloud of small gray points represents instant realizations of the densities numerically found for the exact nonlocal equations. The separatrix connecting this point and the saddle point divides plane {ρ1,ρ2} into attracting basins for the two stable directly seen fixed points. Studying the behavior of the trajectories that start near the separatrix, we can predict several results of numerical modeling of the complete equations and, in the final analysis, the properties of real systems. In particular, we can easily predict the role of the source of noise. If the noise is strong, the phase trajectories can pass both above and below the separatrix, irrespective of the scenario according to which the front arrives at the given point. Within a certain time interval after the arrival of the front, both densities ρ1,2 increase essentially simultaneously and very fast, to which the first maximum in the evolution rate Methods of the Pattern Formation in Numerical Modeling of Biological Problems 231 2 2 2 1 2 ( ) [( / ) ( / ) ]W t t t       (24) corresponds, as depicted in Fig. 4. Near the saddle point, there can be no further increase in ρ1,2. The rate W(t) rapidly decreases. At the same time, there is phase separation in the system, with one of the densities, ρ1 or ρ2, expelled from the given region in space. This is followed by a sharp increase in W(t) accompanied by a rapid buildup of the remaining component, a process that is slowed down near one of the stable fixed states. Evolution rate W(t) as a function of time is shown in Fig. 6. Fig. 6 Evolution rate W(t) as a function of time. Good correlation with the phase portrait (shown in the insert) is clearly seen. The particular trajectory for which W(t) is calculated is marked by red color. Low evolution rate corresponds to the vicinities of the fixed points and “large river” leading to final static configuration. Discrete set of black points represent the mean rate numerically found for the exact nonlocal equations. The characteristic double-humped curves representing the evolution rate are indeed observed, when the complete system of equations is solved numerically. In accordance with the physics of the problem, the initial condition is selected in the form of a narrow strip of density ρ1(r,t), near one of the boundaries of the two-dimensional system. Here the process of generation and separation of the densities ρ1,2(r,t) is accompanied by the formation of characteristic dendritic spatial distributions of both densities. The two densities are generated simultaneously in the vicinity of the front. However, in the absence of noise, the initial density distribution leads to the formation behind the front of a completely filled region. Numerical solution shows that further with time the front becomes discontinuous, and the expanding ordered region is transformed into a fractal. In Fig. 5, we use a copper scale to show a characteristic fragment of the system with a distribution of the total density that emerges at the intermediate stage in the transition from the homogeneous growth to the fractal growth. Contaminated regions and 232 A.E. FILIPPOV, S.N. GORB regions of active growth are clearly visible. These are characterized by intermediate values of the densities which correspond to sections in Fig. 5 with intermediate shades of copper color. Fig. 7 Two intermediate stages of the porous surface layer formation (from subplot (a) to (b)). With the time, such system can develop very complex structure with not simply fractal surface, but also with lacunas of a negative curvature, which effect in a strong increase of adhesion of the soft tissue contacting such a surface. In comparison to the study of growth processes with models with only one fluctuating variable, the novel property of the system lies in the possibility of the emergence of voids behind the front, i.e. regions filled with neither of the two components. When there is only one field, i.e. ρ1(r,t) (such lacunae must be thought of as being filled with the other, „„contaminating‟‟ component of density, ρ2(r,t), which in this case is not explicitly present in the equations. The model contains additional information about the second field ρ1(r,t), which makes it possible to distinguish between the regions occupied by ρ2(r,t) and the voids. The mechanism of void formation is clearly seen in Fig. 7a, which for a small fragment of the front depicts a typical growth sequence. Three characteristic moments in time are singled out: the emergence of a dense initial layer, the emergence of the first dendritic protuberances, and the collapse of the first internal pores in a structure with a total density ρ(r,t)=ρ1(r,t)+ρ2(r,t). Well-formed voids are clearly visible in Fig. 7b. Void formation is closely related to the ability of the „„contaminant‟‟ ρ2(r,t) to block the active sections of the front ρ1(r,t) and, at least in principle, to terminate the growth process. As the pores collapse, the front usually continues to move in both directions. Naturally, the outer boundary ρ1(r,t) is almost insensitive to the presence of a pore blocked somewhere inside the system and continues its forward motion. The inner boundary ρ1(r,t) surrounding a pore is qualitatively similar to the outer boundary and can move „„back,‟‟ up to the point, where it is completely blocked by sections with ρ2(r,t). The scenario of the evolution of the system turns out to be exceptionally multifaceted and Methods of the Pattern Formation in Numerical Modeling of Biological Problems 233 under a slight variation of the coefficients of the model makes possible a reproduction of very realistic configurations of the densities ρ1(r,t) and ρ2(r,t). 6. FORMAL GENERATION OF FRACTAL SURFACES It was an example of creation of the fractal surface by growth. Now, let us turn to more formal methods of numerical generation of surfaces. Depending on a particular problem, we use at least two ways to generate rough surface that is close to the fractal one. One is to apply Fourier expansion with the harmonics restricted between maximal and minimal wavelengths. Depending on the number of the waves and both limits, the surface will be more or less close to the real fractal one and include or not irregularities with small sizes [10]. Another way is to use random deposition of the Gaussians. Depending on their number, height dispersion, minimal and maximal widths, the surface will be closer to fractal or not. Many of biological surfaces have a structure made as a pattern of close objects with some specific scale. They can be more conveniently simulated by a set of deposed objects. That is why, for our studies we usually choose the second method, which easily provides such possibility. In many cases, such surface appears as a result of random deposition of some particles, balls, block of some other objects. We definitely admit that the formal model completely ignores the fine structure of adhesion-like weakening at motion, fracture and restoration of the bonds, so on. Often just only two basic mechanisms are involved into our minimalistic models: attraction to the surface and energy dissipation. The first one actually restricts the scales, which we include to the surface by the irregularities close to the “size of contact points” (radius of attraction in the model). In particular, this limitation corresponds to a situation with adhesive terminal elements (spatulae) of attachment organs of different animals, which were studied elsewhere [11]. The spatulae practically ignore irregularities much smaller or larger than their size. As result, such surfaces can be treated as almost flat. That is why we normally use here the surfaces with the irregularities comparable to the spatula. Such deposition with fixed interval of scales can be done either naturally or artificially. We will not concentrate on any specific case and discuss only a couple of more or less abstract examples. Despite of formal mathematical generality such approach helps in developing a simple, but realistic numerical model that might be useful, for example, in the study of animal spatula interaction with various substrate profiles. Formal models can well explain results of experimental studies, and also predict adhesion of animals to the real surface profiles depending on the dimension and stiffness of the spatula. Previously, we studied such biological systems using numerically generated surfaces and compared the results obtained with the real 3D surface profiles obtained experimentally by depositing of spheres. In the purely numerical approach, different roughness of the surfaces can be achieved by depositing of spheres with defined radius on originally flat surface. To simulate this process numerically, we organized our procedure as follows. At every given radius R , we have taken an array of equal spheres, or circles in two- dimensional (D=1+1) case and placed them successively into the randomly chosen positions xn, where n=1,2,…,N. The number of the spheres N=||L/R|| in the array was taken to be an integer value corresponding to the total number of the circles of radius R 234 A.E. FILIPPOV, S.N. GORB which are necessary to cover the system of the length L uniformly. Each sphere was added virtually to the corresponding segment of the surface: 2 2 ( ; ) ( ) n n n y x x R x x   . (25) Due to the random deposition, the segments can fall either on the top of already existing coverage or onto the empty space. As a result of such a deposition, the total surface gradually accumulates all the segments 1 ( ) ( ; ) N n n n y x y x x    , which are generally speaking placed one on top of another one. This procedure can be easily generalized for D=2+1 surfaces. Typical realizations of such surfaces for different radii R are presented in the Fig. 8. Fig. 8 Consequent stages of the formation of rough surface by the random deposition of the hard balls with fixed radius on the originally flat surface In many cases specific form of the deposed particles is not very important (or maybe not important at all). In such a case, one can treat random deposition procedure as an absolutely abstract way of constructing desirable surface. One of the simplest ways to do this in this case is to model rough surface Z(x,y) by a random deposition of the Gaussians: 2 2 2 ( , ) ( , ,{ , }) exp[ (( ) ( ) ) / ] n n n n n n n n n Z x y G x y x y a x x y y w       . (26) The advantage of this approach is its mathematical transparency. One can regulate everything by simple manipulation with varied positions, amplitudes and widths of the Gaussian functions. The only limitation is generally that the characteristic scale of the artificial surface irregularities should be adjusted to be smaller or comparable with the scale of the elements interacting with the surface (spatulae of animals, for example). In particular, the typical distance between the hills and valleys of the randomly accumulated surface Z(x,y)=Gn can be regulated by the number of Gaussians. An additional convenience of the method is in the fact that one can even not control the amplitude of the asperities during accumulation of the sum Gn. The amplitude of roughness after accumulation is finally regulated by the normalization Z(x,y) → A(Z(x,y)−min(Z))/(max(Z)−min(Z)), where desirable amplitude A can be chosen from the limit of the flat surface A=0 to the values comparable with the characteristic lengths of the system under consideration. Different variants of the rough surface generated by the random deposition of Gaussians are shown in Fig. 9. Methods of the Pattern Formation in Numerical Modeling of Biological Problems 235 Fig. 9 Different variants of the rough surface generated by the random deposition of Gaussians. The number of Gaussians used to generate the surface grows from left to right. One of the widely known reasons to generate numerical surface or potential U is its application to tribology. Let us mention, in this context, a commonly used Prandtl- Tomlinson model that has proven to be successful in describing shear response in surface force apparatus configuration. This model describes the lateral motion of a driven plate and has the following form: 2 2 / / ( ) / ( ) 0x t x t U x x K x Vt          (27) Here a driven plate of mass M and the center-of-mass coordinate x is pulled by a spring of a spring constant K. In standard Tomlinson equation that is normally applied to study atomically flat surfaces, U(x) is the effective periodic potential U(x)=U0cos(2x/b) experienced by the plate due to the presence of an embedded system. In some sense, such simple potential is also particular realization of the surface pattern appearing on microscopic level, where its creation is regulated by the process of crystallization, which normally leads to an ideally periodic crystal lattice. Below for general theoretical study, we use dimensionless units normalized to characteristic microscopic scales and energies of the particular system under consideration. Parameter , as usually, is responsible for the dissipation. The spring is connected to a stage which moves with a velocity V. Tomlinson equation was well supported microscopically, but usually people do not extend its application even to the mesoscopic scales. However, as we saw in many examples of the formation of many surfaces (may be better to say all the surfaces which we meet in normal life), the surface is dictated by frozen kinetics, which leads to the relief, which is quite far from an ideal periodic structure. As a result, one partially knows how to relate the parameters that appear in Eq. (27) to the microscopic characteristics of the system, but can not correlate them with the macroscopic friction. This problem is not solely limited to the friction studies, but appears also when dealing with adhesion. In the previous papers, we have seen how important it is for the study of attachment devices of animals, which have to adapt to real surfaces also on intermediate mesoscopic scales, but not only at the microscopic one, where pattern either can be approximately treated as a flat or even simply periodic [12]. So, what becomes a major problem here is that the mesoscopic structure of frictional surfaces is of a fractal character and thereby cannot be characterized by a certain wave vector (or even few ones), like it is normally used in applications of the Tomlinson model. On the other hand, it could be inconvenient or simply senseless to reproduce each time one of the physical or chemical processes described here. It is important to know, 236 A.E. FILIPPOV, S.N. GORB how such surfaces appear, or maybe it becomes important, if the studied surface really has specific substructure with the periodicity important for the modeling. However, if we can forget about underlining process and accept a naive hypothesis that the majority of the surfaces are simply fractal, it is natural to extend the model into a directly generated pattern according to very formal mathematical definition of the fractal. Let us consider a fractal potential of the form: 2 1 0 ( ) ( ) cos( ) q q U x U dqc q qx   , ( )c q q    (28) Here q1 and q2 are characteristic cut off wave vectors and δ(x) is the random phase that we assume δ-correlated <δ(q)δ(q‟)>=2πδ(q-q‟). For the majority of physically interesting systems, index β is close to β≈0.9. Below we will keep this value for definiteness. For further study, it is convenient to go over to a discrete representation of the integral in Eq. (2) ∫dqc(q)→Σ with a discrete step between the wave vectors Δq determined by the smallest vector q1 corresponding to an inverse maximal length lmax of the system, which equals normally to its size lmax=L . Total number Ntot of the terms in the sum is given by Ntot=q2/q1≡q2/Δq. The discrete approach is natural for our further numerical studies. It allows us to adjust the potential to different scales by a number of the modes included into summation Nmodes