Piotr Dziekan 
  
prof. nadzw. IChF PAN Bogdan Nowakowski 
  
Annie Lemarchand 
  
Dynamics of far-from-equilibrium chemical systems: microscopic and mesoscopic approaches

Many nonlinear systems under non-equilibrium conditions are highly sensitive to internal fluctuations. In this dissertation, stochastic effects in some generic reaction-diffusion models are studied using two approaches of different precision. In the mesoscopic approach, evolution of the system is governed by the master equation, which can be solved numerically or used to set up kinetic Monte Carlo simulations. On the microscopic level, particle computer simulations are used. These two stochastic approaches are compared with deterministic, macroscopic reaction-diffusion equations.

In the Introduction, key information about the different approaches is presented, together with basics of nonlinear systems and a presentation of numerical algorithms used.

The first part of the Results chapter is devoted to studies on reaction-induced perturbation of particle velocity distributions in models of bistability and wave front propagation. A master equation including this perturbation is presented and compared with microscopic simulations.

The second part of the Results deals with pattern formation in reaction-diffusion systems in the context of developmental biology. A method for simulating Turing patterns at the microscopic level using the direct simulation Monte Carlo algorithm is developed.

Then, experiments consisting of perturbing segmentation of vertebrate embryo's body axis are explained using the Turing mechanism. Finally, a different possible mechanism of body axis segmentation, the "clock and wavefront" model, is formulated as a reaction-diffusion model.

I. Introduction

Every closed physical system evolves towards the equilibrium state, in which all thermodynamic flows vanish and, consequently, the system is highly homogeneous [1]. Yet everywhere around us we observe stable spatial structures, and even evolution towards higher complexity, which seems to be in contrast to this principle. To understand this kind of self-organization of complex systems, we have to shift our attention from the vicinity of equilibrium to systems that are far away from it [2]. To be stable, such an out-of-equilibrium state has to be sustained by flows of matter or energy between the system and its surrounding [3,4]. Almost every interesting complex system seen in Nature, from living cells to the Earth as a whole, can function only thanks to such flows.

Removing a system from equilibrium results in appearance of some kind of order in it, which is related to the decrease of its entropy [5,6]. Qualitatively new behaviour and structures can appear only in nonlinear systems, that is systems in which thermodynamic flows are not linear functions of thermodynamic forces. Chemical reactions very often have this nonlinear property. This work deals with chemical systems, but the methods used can be easily applied to other dynamical systems, like population dynamics or disease spreading.

Real chemical systems are characterised by very high level of complexity, usually having many types of species and interactions between them. This causes two major problems in studying them. One of them is the difficulty in experimentally determining all necessary parameters, like reaction rate constants and diffusivities. The other prob-lem lies in understanding dynamics of such a complicated network of interactions. For these reasons simplified models are usually used. They capture crucial features of real systems and at the same time allow detailed analysis. This work is focused on reactiondiffusion models, i.e. models in which only chemical reactions and diffusion of particles are included. Every model used is based on a scheme of elementary reactions.

The basic method of description of reaction-diffusion systems is by deterministic macroscopic equations. Often these equations give an oversimplified view of the phenomenon, therefore in this work special attention is given to other, more realistic methods. On the mesoscopic level, dynamics is governed by the master equation, which is intrinsically stochastic (Sec. 1.2). The third method used are the microscopic computer simulations, described in detail in Sec. 4.

Subject of my research can be divided into two general categories. Section 5 deals with the influence of perturbation of particle velocity distributions on the dynamics of chemical systems. The second subject is formation of one-dimensional periodic patterns in the context of biology (cf. Sec. 6).

Different description levels of dynamics of chemical systems

This work deals with dynamics of chemical reaction-diffusion systems, i.e. temporal evolution of a mixture, in which chemical reactions and diffusive transport take place. An elementary reaction can be either caused by internal processes within a single molecule, when we speak of a first-order reaction, or a result of collision with another molecule, what we call a reaction of second order. An instantaneous collision of more than two molecules is extremely unlikely, therefore reactions of higher order are treated as series of successive binary collisions taking place in a short time interval [7]. Rigorously, state of a chemical mixture should be described by a discrete number of particles of each type. On the other hand, it is more convenient to work with species concentrations, i.e. number of particles per unit volume, which are continuous variables. Such approximation is valid in the case of dense systems and high molecular populations. In addition, when number of molecules is large, we can assume that species concentrations change in a deterministic fashion. Description based on these two assumptions, i.e. continuous variables and deterministic dynamics, will be called a macroscopic one.

The second, more accurate level used is the mesoscopic description. It is based on the master equation and takes into account the discreetness of molecular populations as well as the stochastic nature of chemical reactions. It is worth to note that there are stochastic equations for continuous variables, like the Fokker-Planck equation, but their derivation requires some additional assumptions about the properties of fluctuations and they may give wrong results in systems with complex dynamics [8][9][10][11]. Using the mesoscopic approach reveals new, qualitatively different behaviours of the system and is necessary when analysing small systems, like living cells, in which fluctuations play a significant role.

The description method that is closest to reality relies on tracking positions and momenta of every single molecule, what we call the microscopic level. Due to an extremely large number of variables, no analytic solutions are known on this level and the only tool available are numerical simulations. Algorithms of such simulations are described in Sec. 4.3.

Macroscopic approach

In the macroscopic approach, dynamics of a system containing N reactants is governed by a system of N differential equations in the form [12]:

∂ t c i (x,t) = f i (c 1 ,c 2 ,...,c N ) -∇ • j i , (I.1)
where c i (x,t) is the concentration of the i-th reactant at position x and at time t, j i is the diffusive flux and f i (c 1 ,c 2 ,...,c N ) = j a ij g j (c 1 ,c 2 ,...,c N ) is the total rate of change of concentration of species i due to chemical reactions [13,14]. The coefficient a ij is equal to the change in the number of type i particles caused by a single j-th reaction.

In the case of elementary reactions, i.e. reactions that do not include any intermediate products, reaction rate g j is given by the following equation [15,16]:

g j = k j c s 1j 1 c s 2j 2 ...c s N j N , (I.2)
where k j is the reaction rate constant and s ij the stoichiometric coefficient associated with i-th species in the scheme of j-th reaction. For example, in case of reaction A +

B k 1 --⇀ ↽ -- k -1 2A, we have f A = -f B = k 1 c A c B -k -1 c 2 A .
In general, diffusive flux of i-th species j i depends on spatial inhomogeneities in concentrations of all species of the mixture [17]. However, in many cases cross-diffusion is negligible and we can assume that diffusive flux of species is proportional only to the gradient of its concentration. This approximation is expressed by the Fick's law of diffusion j i = -D i ∇c i , where D i is the diffusion coefficient of species i. This linear dependence of the diffusive flux on the gradient of concentration is correct if Knudsen number, which is the ratio of the mean free path of molecules to a characteristic length of the inhomogeneities of concentrations, is less than 0.1 [18]. If the diffusion coefficient is constant, equation (I.1) takes the form:

∂ t c i (x,t) = f i (c 1 ,c 2 ,...,c N ) + D i ∆c i . (I.3)
Except for some simple models, systems of equations of the type (I.3) have no analytic solutions. To obtain a result, numerical methods have to be employed, which will be discussed in Sec. 4.1.

Mesoscopic approach -master equation

Natural variables for a chemical system are species populations (number of particles of each type), which we denote N i (t). For now, let us focus on homogeneous systems. The state of a system is given by the vector N , which consists of populations of all species.

Since on this description level we want to include intrinsic fluctuations, we are interested in the temporal evolution of probability distribution of system states P (N ,t). Reactions cause discrete changes of species populations and dynamics of a chemical system can be seen as a birth and death process. Assuming that probabilities of subsequent reactions are independent, i.e. reaction is a Markov process, an equation can be derived for the time derivative d t P (N ,t), which is called the master equation [5,19,20]:

d t P (N ,t) = j
w j (Nr j )P (Nr j ,t)j w j (N )P (N ,t).

(I.4)

The vector r j represents the change of the vector of state due to a single j-th reaction.

By w j (N ) we denote the probability of transition per unit time from the state N to the state N + r j caused by reaction j. The first sum in the master equation represents flux of probability towards the state N , while the second sum is the total flux of probability from this to other states. These sums are taken over all the reactions.

Transition probabilities for elementary second-order reactions are proportional to the number of pairs of particles that can collide [5]. For example, in the model A+B

k 1 --⇀ ↽ -- k -1
2A, already used in Sec. 1.1, they are:

w 1 (N A ,N B ) = k m 1 N A N B (I.5) w -1 (N A ) = k m -1 2 N A (N A -1). (I.6)
The transition probabilities have similar form to the rate equations (I.2), but the constants k m 1 and k m -1 are different from the macroscopic rate constants, because they appear in equations for populations instead of concentrations. In addition, there is a combinatorial factor 1 2 in the expression for w -1 . We have k m 1 = k 1 /V and k m -1 = 2k -1 /V , where V is the system volume. The master equation is more fundamental than the macroscopic reaction-diffusion equation, which can be derived from it [19].

The deterministic number of particles of species i is defined as the mean value of N i : N i = N N i P (N ). If we write the master equation in its integral form, we can obtain the macroscopic equation as the first approximation of d t N i .

The equation (I.4) was defined for a homogeneous system, so it includes fluctuations only of the total number of particles of each species. Often, by extending results for systems close to equilibrium, it is assumed that amplitude of the internal fluctuations is proportional to the reciprocal of the square root of the number of particles. Therefore probability of large fluctuations that would cause qualitative changes in the whole system is very low for systems of macroscopic sizes. Such a global approach to fluctuations is valid for linear systems. In nonlinear systems however, it is no longer applicable, because in these systems local fluctuations may cause qualitative changes in some small part of the system [5]. Next, this change could spread on the rest of the system due to diffusion. Probability of large, local fluctuations is much greater than that of a large global fluctuation, because of the smaller number of particles involved. Therefore it is necessary to use a master equation that would include spatial inhomogeneities in the system. In addition, such equation would be useful in the case of non-homogeneous initial conditions.

The system is divided into spatial cells, between which molecules can diffuse. For simplicity, let us consider a system containing only one type of particles. We are interested in the probability distribution P ({N }), where {N } is a set of populations of these molecules in all spatial cells. The cells have to be small enough to be considered as homogeneous despite possible macroscopic inhomogeneities. On the other hand they cannot be too small, because transport between them is described by diffusion and not free motion. A master equation for inhomogeneous systems that uses discretization of space has the form:

d t P ({N }) = k j
[w j (N k -r j )P (N k -r j ) -w j (N k )P ] The first term of Eq. (I.7) is the probability flux due to chemical reactions, as in the homogeneous case, and the second one represents the probability flux due to diffusion. Form of the transition probabilities associated with diffusion depends on space discretization. In general, transition probabilities can be derived from Fick's laws of diffusion using the finite difference method [21]. In a one-dimensional system and with cells of the same size, divergence of gradient of population in the k-th cell can be calculated as 2 , where ∆x is the cell size [5]. The respective transition probability for diffusion is

(N k-1 + N k+1 -2N k )/(∆x)
d k→l (N k ) = DN k ∆x 2 δ k±1,l , (I.8)
where δ i,j is the Kronecker delta.

The master equation is an intermediate approach between the reaction-diffusion equations and the microscopic description. Except for the simplest situations, it does not have an analytic solution. Still, numerical simulations of system dynamics based on the master equation are possible using the kinetic Monte Carlo (KMC) algorithm (cf.

Sec. 4.2).

Nonlinear systems

The aim of this Section is to present a brief introduction to the dynamics of nonlinear systems [22][23][24]. The diversity of qualitatively different behaviours of a dynamical system increases with the number of variables describing it. Work presented in this thesis is devoted to effects, that are already present in one-variable (multistability and wave fronts) or two-variable (chemical oscillations, Turing patterns) systems. Therefore systems with higher number of variables, which can show chaotic behaviour [25], will not be presented.

In chemical systems the variables are usually concentrations of species. At the macroscopic level, evolution of the system is described by a set of deterministic equations for the time derivatives of concentrations that have the form (I.1). Concentrations values for which all time derivatives are equal to zero, are called stationary states. Except some critical cases, which do not play an important role in real systems, a stationary state can be either stable, that is when every small perturbation of it vanishes, or unstable, when a small perturbation causes a transition to a different part of the phase space (space of all variables of the system). In the deterministic approach, without an external perturbation, the system is unable to leave a stationary state, stable or not.

In the mesoscopic approach, the system is characterised by a stationary probability distribution, for which d t P ({N }) = 0. Irrespective of the initial probability distribution, in time it will develop towards the stationary one. The stationary distribution has local maxima for species populations corresponding to the stable stationary states of the macroscopic equation, but the instantaneous state of the system constantly fluctuates around this maxima. Local minima of the stationary probability distribution correspond to unstable stationary states.

In general, a point or a set of points in the phase space, towards which a dynamical system converges over time, is called an attractor. Regardless of the initial state, system will reach one of the attractors for t → ∞. Part of the phase space from which the system evolves towards a given attractor is called a basin of attraction of this attractor.

In a two-dimensional phase space, if there is more than one attractor, their basins of attraction are separated by separatrices, which may cross through unstable stationary states, or unstable limit cycles (limit cycles will be defined in Sec. 2.1.2). Line in the phase space for which one of Eqs. (I.1) is equal to zero is called a nullcline. The stationary states are located at the crossing of all nullclines. Stability of a stationary state can be determined by the linear stability analysis. This method, along with most important types of stationary states, will be presented in Sec. 2.1.

Homogeneous systems

Homogeneous systems are systems in which spatial differences in species concentrations are negligible. An experimental example of a homogeneous system are reactions performed in stirred reactors. In that case, concentrations depend only on time and transport effects are not included in the equations. State of a mixture of N species is given by the vector c = [c 1 ,c 2 ,...,c N ], where c i is the concentration of i-th species. At a macroscopic level, dynamics is governed by a set of N differential equations:

d t c 1 = f 1 (c) . . . (I.9) d t c N = f N (c).
Stationary states are found by solving d t c = 0. Stability of a stationary state c s can be determined by expansion of Eqs. (I.9) into a Taylor series around that state [23,24].

Since we are interested in stability with respect to small perturbations u = cc s , we only keep linear terms of the expansion. Including f i (c s ) = 0, we obtain: [26,27]:

d t u T = J s u T , J =       ∂f 1 ∂c 1 . . . ∂f 1 ∂c N . . . . . . . . . ∂f N ∂c 1 . . . ∂f N ∂c N       , (I.
• λ 1 ∈ R, λ 2 ∈ R, λ 1 = λ 2 , λ 1 < 0, λ 2 < 0
Stable node -the system will exponentially in time reach the stationary state along one of the eigenvectors. Figure I.1(a) shows trajectories for the case λ 2 < λ 1 , in which the system close to the stationary state moves along the eigenvector v 1 .

•

λ 1 ∈ R, λ 2 ∈ R, λ 1 = λ 2 , λ 1 > 0, λ 2 > 0
Unstable node -the system will exponentially in time leave the stationary state along one of the eigenvectors.

• λ 1 ∈ R, λ 2 ∈ R, λ 1 = λ 2 , λ 1 < 0, λ 2 > 0
Saddle point -an unstable stationary state shown in Fig. I.1(b). Only small perturbations exactly along the eigenvector v 1 will vanish. Separatrices of the system pass through a saddle point along eigenvectors.

• λ 1 ∈ C, λ 2 ∈ C, λ 1 , λ 2 = α ± iβ, α < 0
Stable focus -solution of the set of equations (I.10) includes terms proportional to exp[±iβt], so the system will reach the focus in an oscillatory fashion, as shown in •

λ 1 ∈ C, λ 2 ∈ C, λ 1 , λ 2 = α ± iβ, α > 0
Unstable focus -system will leave the stationary state in an oscillatory fashion.

• λ 1 ∈ C, λ 2 ∈ C, λ 1 , λ 2 = ±iβ
Centre -elliptic trajectories around the stationary state (cf. A simple model of a bistable system (having two stable states) is the Schlögl model [28,29]:

2X + A k 1 --⇀ ↽ -- k -1 3X (I.11) X k 2 --⇀ ↽ -- k -2
B, (I.12)

where concentrations of A and B species are constant. This assumption is valid if concentrations of A and B are much larger than that of X, or if they are maintained constant through coupling to external reservoirs, like in a continuous stirred tank reactor.

System dynamics is governed by a single polynomial of third degree:

d t X = -k -1 X 3 + k 1 AX 2 -k 2 X + k -2 B, (I.13)
where species concentrations are in italics. This equation can have up to three roots X 0 , X 1 , X 2 , which are two stable stationary states with an unstable one between them.

In a deterministic approach the system can remain in one of the stable states and transitions between them are not possible. Additional effects can appear if stochasticity is included. The stationary probability distribution has two local maxima N X0 and N X2 , corresponding to the stable states. The unstable stationary state corresponds to a local minimum N X1 . The stationary probability distribution satisfies a condition that there are no net probability flows between states. If transitions only between states with populations differing by 1 are possible, like in the case of the model (I.11,I.12), and if an arbitrary value of P (0) is assumed, probabilities of the other states can be found from an iterative formula w + (N X )P (N X ) = w -(N X + 1)P (N X + 1), where w + and w -are probabilities of transition to a state with population greater and less by 1, respectively.

The system may remain for a long time in vicinity of one of the maxima, but a significantly large fluctuation may cause a transition to the other maximum. The mean time of such passage can be calculated using transition probabilities from the master equation (I.4) [19]:

τ 0→2 = N X2 -1 ν=N X0 ν µ=0 1 w + (µ) ν ξ=µ+1 w -(ξ) w + (ξ) (I.14) τ 2→0 = N X2 -1 ν=N X0 µ≥ν 1 w -(µ) µ ξ=ν+1 w + (ξ) w -(ξ) , (I.15)
where τ 0→2 is the mean time of passage from the state N X0 to the state N X2 and τ 2→0 is the mean passage time in the reverse direction. If the product at the end of Eqs. (I.14,I.15) does not contain any terms, it is considered to be equal to one.

An example stationary probability distribution P (N X ) of a bistable system is presented in Fig. I.3. The closer the stable states are to each other and the higher is probability of the unstable state, the shorter is the mean passage time. 

Chemical oscillations

In systems of at least two variables a new type of attractor may appear, called a stable limit cycle. It is a closed path in phase space along which the system travels and which attracts surrounding trajectories. In time domain it is represented by regular, sustained oscillations of species concentrations. General conditions for the existence and number of limit cycles are not known. However, the Poincaré-Bendixson theorem states that every simply connected subspace of a two-dimensional phase space that does not contain stable states and to which surrounding trajectories enter, contains a limit cycle [30].

One of the basic ways how a stable limit cycle can appear is through the supercritical Hopf bifurcation [31]. It takes place when, due to a change in parameters, real part of the eigenvalue of the Jacobian changes its sign from negative to positive, while the imaginary part of it remains nonzero. Therefore it is a transition from a stable to an unstable focus.

Let us consider a two-variable system:

d t c 1 = f (c 1 ,c 2 ), d t c 2 = g(c 1 ,c 2 ). (I.16)
We assume that it has only one stationary state (c s 1 ,c s 2 ). To simplify analysis of its stability we introduce new variables

x = c 1 -c s 1 , y = c 2 -c s 2 .
In new variables the stationary state is at the axis origin x = 0, y = 0. Linearizing Eqs. (I. 16) around this state we get:

d t x = f x | (0,0) x + f y | (0,0) y (I.17) d t y = g x | (0,0) x + g y | (0,0) y. (I.18)
Stability of the stationary state depends on eigenvalues of Jacobian of the above set of equations. They can be calculated from the characteristic polynomial of the Jacobian in the stationary state:

det (J s -λI) = λ 2 -λtr(J s ) + det(J s ) = 0. (I.19)
Eigenvalues, which are roots of the characteristic equation, are:

λ 1,2 = tr(J s ) ± tr 2 (J s ) -4 det(J s ) 2 . (I.20)
Therefore the state (0,0) is an unstable focus if: tr(J s ) > 0, tr 2 (J s ) < 4 det(J s ).

(I.21)

A supercritical Hopf bifurcation occurs under such a change of parameter values, that trace of the Jacobian becomes positive and the condition det(J s ) > 0 is fulfilled. The stationary state becomes unstable and the system will leave it in an oscillatory manner.

Since there are no other stationary states, and a realistic reaction scheme cannot allow concentrations to grow to infinity, the system will reach some other type of attractor.

Here, this attractor will be a stable limit cycle surrounding the unstable stationary state.

It is characteristic for this kind of bifurcation that amplitude of the limit cycle grows as the parameters move away from the bifurcation point. A schematic representation of the transition from a stable focus to a stable limit cycle is shown in The Brusselator is a classical model of oscillations in chemical systems [32]. It includes two species of variable concentrations, X and Y:

A k 1 -→ X (I.22) B + X k 2 -→ Y + D (I.23) 2X + Y k 3 -→ 3X (I.24) X k 4 -→ E. (I.25)
It is a model of reactions in an open system in which species A and B are turned into E and D, respectively. It is convenient to introduce scaled variables:

X ′ = k 3 k 4 X, Y ′ = k 3 k 4 Y, t ′ = k 4 t (I.26) A ′ = k 2 1 k 3 k 3 4 A, B ′ = k 2 k 4 B. (I.27)
In new variables, dropping the primes, kinetic equations become:

d t X = X 2 Y -(B + 1)X + A (I.28) d t Y = -X 2 Y + BX. (I.29)
This model has one stationary state X 0 = A, Y 0 = B A . Jacobian in this state is:

J s =    B -1 A 2 -B -A 2    . (I.30)
From Eqs. (I.21) we know that a supercritical Hopf bifurcation happens for:

B = A 2 + 1. (I.31)
Stable limit cycle is a structurally stable attractor. It means that any perturbation from it will disappear in time and in that manner it is different from oscillations in conservative systems, trajectories of which depend on initial conditions [26,33,34]. Therefore a limit cycle is stable to fluctuations, except fluctuations that are parallel to its trajectory, which will not disappear. This effect is called phase diffusion. Let us consider a large number of identical systems, which initially all are at the same point along limit cycle. With time, phases of their oscillations will drift away from each other and for t → ∞ they will not be correlated at all. This effect is pictured in Fig. I.5. Already after a single period of oscillations states of independent oscillators are dispersed over the whole limit cycle. Moreover it can be seen that the system travels fast in parts of phase space far away from nullclines and it slows down close to them.

Reaction-diffusion systems

In many real systems spatial differences of concentrations may appear. Then the equations that govern dynamics of the system have to include diffusive transport of species. Evolution of a system is given by a set of equations of the type:

∂ t c i (x,t) = f i (c 1 ,c 2 ,...,c N ) + D i ∆c i , (I.32)
where D i is the diffusion coefficient of i-th species.

Wave fronts

Let us consider an inhomogeneous chemical system with dynamics that allows it to remain for a long period of time in one of two different stationary states. We assume that part of the system is in one state, and the rest in the other state. Then a diffusive transport of species will appear at the interface between these areas. This transport, together with chemical reactions, will cause an increase of one of the areas at the expense of the other. The propagating gradient of concentrations is called a wave front. In the simplest case of a one-dimensional system of one variable the equation for dynamics is:

∂ t c(x,t) = f (c) + D∂ xx c. (I.33)
We are looking for a solution that would represent a wave front propagating with constant velocity v in the positive direction of the x axis [35]. Also, we assume that this front does not change shape in time. Such a solution is in the form:

c(x,t) = c(x -vt) = c(z), z = x -vt. (I.34)
After substituting it into Eq. (I.33) we get a second-order ordinary differential equation:

Dd zz c + vd z c + f (c) = 0. (I.35)
The method of solving this equation depends on the form of function f (c). An important case are systems that have one stable stationary state c = 1 and one quasi-stable state c = 0. A quasi-stable state is a stationary state that, from the point of view of dynamics is unstable, but the system cannot leave it through fluctuations, because there are no particles of the proper type. Destabilisation of a quasi-stable state can happen only if some particles are introduced from the outside, for example through diffusion.

The front will propagate from the area with c = 1 over the area c = 0, which imposes following constraints on the function c(z):

lim z→∞ c(z) = 0, lim z→-∞ c(z) = 1. (I.36)
Fronts propagating between a quasi-stable and a stable state can be of two types: pulled or pushed [36]. Pulled fronts are characterised by that the derivative f ′ (c) has maximal value at the stationary state c = 0, which is ahead of the front. Then, front speed is solely determined by dynamics of reactions at its leading edge, where f (c) ≈ 0, therefore we can use a linear approximation of the function f (c), which yields:

Dd zz c + vd z c + f ′ (0)c = 0. (I.37)
This equation has a realistic, non-negative solution if:

v ≥ v c = 2 Df ′ (0). (I.38)
Therefore speed of stable pulled fronts has to be greater than a critical value v c . Width of a front is defined as the absolute value of the reciprocal of its gradient at the inflection point. Assuming that the inflection point is at c = 1/2, front width is w = v/f (0.5), so it is proportional to the front speed. Since the front speed is not strictly defined, it depends on the initial conditions: speed of fronts narrower than v c /f (0.5) will develop to v c , while fronts that are wider will propagate adequately faster [37][38][39].

An important example of a pulled front is the Fisher-Kolmogorov-Petrowsky-Piskunow model (FKPP) [40,41]. It is a general model used in studies of many effects, like propagation of a favourable gene in population or disease spreading. In the context of chemistry it represents propagation of a wave front due to the reaction

A + B k - → 2A
with conservation of the total number of particles: c A + c B = 1. The system is described by one variable c = c A , dynamics of which is:

f (c) = kc(1 -c). (I.39)
The critical speed of the FKPP wave front is v = 2 √ kD and its respective width If the derivative f ′ (c) has a maximum for c > 0, propagation speed is determined by dynamics of reactions not at its leading edge, but inside of the wave front. Such fronts are called pushed fronts. There is no general, analytic expressions for their speed, but it is also bounded from below by a critical value v p , with v p > v c .

Many reaction-diffusion models of real systems admit more than one stable state, so wave front propagation between such states is an important question. Such front, called a trigger front, is qualitatively different from the fronts discussed so far, because between the two stable states there is always an unstable one, which separates their basins of attraction [42]. The simplest model of a trigger front is a function f (c) with three roots: 

c 0 < c 1 < c 2 .
Propagation speed of a stable trigger front is strictly defined and equal to 

v = - c 2 c 0 f (c)dc
v = k -1 D/2(c 0 + c 2 -2c 1 ),
where D is the diffusivity [26]. Direction of propagation is determined by the relative probabilities of stable states. For example, for the distribution showed in Fig. I.3, front would propagate from the more probable state N X2 towards the state N X0 .

Turing patterns

In his work about mechanisms of morphogenesis, that is formation of shape during growth of organisms, Turing showed that chemical reactions coupled with diffusion can cause spontaneous formation of stable spatial patterns of species concentrations [26,43].

These patterns can form in systems with at least two species and it is crucial that these species have different mobilities. A broader theory of the so-called dissipative structures, i.e. structures forming in systems that are maintained far from equilibrium by flows, was developed by Prigogine [5,44,45].

Conditions for formation of Turing patterns will be derived for a simple one-dimensional system of two variables u and v.

∂ t u(x,t) = f (u,v) + D u ∂ xx u (I.41) ∂ t v(x,t) = g(u,v) + D v ∂ xx v. (I.42)
We assume that the state (0,0) is stationary, i.e. f (0,0) = g(0,0) = 0. Length scale is chosen to get D u = 1. Ratio of diffusion coefficients is denoted by d = D v /D u . Using linear stability analysis we want to find conditions, under which this state is stable to homogeneous perturbations, but unstable to spatially periodic perturbations.

Linear stability analysis of two-variable homogeneous systems was done in Sec. 2.1.2.

From Eq. (I. 19) we conclude, that the state (0,0) is homogeneously stable if determinant of Jacobian in this state is positive and its trace is negative:

trJ | (0,0) = (f u + g v )| (0,0) < 0, |J | = (f u g v -f v g u )| (0,0) > 0. (I.43)
Now we will perform linear stability analysis of the full system of equations, including diffusion. Let us put the solution in the form of a linear combination of harmonics

u(x,t) = i=∞ i=1 a i exp(λ i t + ik i x) (I.44) v(x,t) = i=∞ i=1 b i exp(λ i t + ik i x). (I.45)
Linearizing functions f and g around point (0,0) and including (I.44,I.45), the set of equations (I.41,I.42) takes the form:

λ i u = f u (0,0)u + f v (0,0)v -k 2 i u (I.46) λ i v = g u (0,0)u + g v (0,0)v -k 2 i dv, (I.47)
for every wavenumber k i . Coefficients λ i are eigenvalues of the matrix

   f u -k 2 i f v g u g v -k 2 i d    (0,0) , (I.48)
therefore they can be found solving:

λ 2 i + λ i k 2 i (1 + d) -(f u (0,0) + g v (0,0)) + h(k 2 i ) = 0, (I.49) h(k 2 i ) = dk 4 i -(f u (0,0)d + g v (0,0))k 2 i + f u (0,0)g v (0,0) -f v (0,0)g u (0,0). (I.50)
For every wavenumber there are two eigenvalues λ ± i . If real part of one of them is positive, the stationary state is unstable to periodic perturbations of respective wavenumber.

As a result, a stable spatial structure will develop. This can happen if coefficient next to the first power of λ i in above equation is negative, or if h(k 2 i ) < 0. The first condition cannot be fulfilled, because from (I.43) we have f u (0,0) + g v (0,0) < 0.

(I.51)

Recalling that, according to the second condition (I.43), we have f u (0,0)g v (0,0)f v (0,0)g u (0,0) > 0, we conclude that the second possibility, i.e. h(k 2 i ) < 0 for some k i > 0, can happen only if the coefficient next to k 2 i is negative: 

f u (0,0)d + g v (0,0) > 0. (I.
f u (0,0) > 0, g v (0,0) < 0. Moreover, the condition (f u g v -f v g u )| (0,0) > 0 (cf.
Eq. (I.43)) imposes that f v (0,0)g u (0,0) < 0, which allows two different choices of signs of these derivatives, that give qualitatively different patterns. If f v (0,0) > 0, g u (0,0) < 0, periodic patterns of the two species have opposite phases, so minimum of concentration of one of them will be at the position of maximum of the other. On the other hand, for f v (0,0) < 0, g u (0,0) > 0 the structures will be in phase [26].

The condition (I.52) is necessary, but insufficient to get h(

k 2 i ) < 0. A sufficient condition is h(k 2 m ) < 0,
where k 2 m is the position of minimum of function h(k 2 i ):

k 2 m = 1 2 f u + g v d , (I.53) 
from which we obtain:

(df u (0,0) + g v (0,0)) 2 4d > f u (0,0)g v (0,0) -f v (0,0)g u (0,0). (I.54)

Knowing rate constants, from the above inequality we can derive a condition for the minimal ratio of diffusivities. We denote this critical ratio by d c . It can be found by solving the equation

d 2 c f 2 u (0,0) + 2(2f v (0,0)g u (0,0) -f u (0,0)g v (0,0))d c + g 2 v (0,0) = 0. (I.55)
For every wavenumber fulfilling h(k 2 i ) < 0 there is at least one eigenvalue with positive real part. A relation between the greater one of Re λ ± i and wavenumber k i is called the dispersion relation.

One of the simplest models of Turing pattern formation is the Gray-Scott model [46,47]. It postulates two reactions taking place in a flow reactor

2X + Y κ 1 -→ 3X (I.56) X κ 2 -→ P. (I.57)
Overall, they account for production of species P from substrate Y . A flow φ is imposed in the reactor, that supplies species Y at rate φY 0 and removes species X, Y and P. The reaction-diffusion equations are

∂ t X = X 2 Y -(φ + κ)X + ∆X (I.58) ∂ t Y = -X 2 Y + φ(Y 0 -Y ) + d∆Y. (I.59)
The time and spatial scales have been selected to obtain κ 1 = 1 and D X = 1, resulting in κ 2 = κ and D Y = d. The system may have up to three stationary states

(0,Y 0 ),(X -,Y -),(X + ,Y + )
, where

X ± = φY 0 ± √ ∆ 2(φ + κ) , Y ± = φ + κ X ± , ∆ = (φY 0 ) 2 -4φ(φ + κ) 2 . (I.60)
Let us consider formation of Turing pattern from the state (X + ,Y + ) in a one-dimensional Gray-Scott system of length L. The Jacobian in this state is

J s =    φ + κ X 2 + -2(φ + κ) -X 2 + -φ    . (I.61)
The choice of boundary conditions is of key importance. We choose Neumann condi-

tions: [48] ∂ x u(0,t) = ∂ x v(0,t) = ∂ x u(L,t) = ∂ x v(L,t) = 0, (I.62)
that limit possible structures to those with wavenumbers equal to k n = n π L , n = 1,2,.... 

k min k max 0.6 0.8 1 1.2 k -0.1 0 0.1 0.2 Re(λ)

Aim of the Thesis

Effects that are characteristic for far-from-equilibrium nonlinear dynamical systems that were presented in Sec. 2 are highly sensitive to internal fluctuations. Propagation speed of pulled wave fronts depends on dynamics of reactions at the leading edge, where species populations are low, so fluctuations have a large impact on the development of the whole wave front. In bistable systems, spontaneous transitions between stable states are due to fluctuations. In the case of chemical oscillators, fluctuations cause phase diffusion. Turing patterns can develop from a spatially homogeneous stationary state only, if a stochastic perturbation moves the system away from it. Another reason why it is important to include fluctuations in our description is that the above-mentioned effects play an important role in living cells [26]. Chemical wave fronts propagate information within and between cells, limit cycles act as biological clocks and cells can remain in one of many stationary states. Finally, Turing patterns account for formation of spatial structures from homogeneous groups of cells. Since species populations within living cells are low, analysing these nonlinear effects in the biological context requires fluctuations to be taken into account.

The aim of this work is to investigate influence of internal fluctuations on behaviour of selected nonlinear systems through analysis on the mesoscopic and microscopic levels. The deterministic, macroscopic approach will be used as a reference. Complexity of chemical systems prohibits solving the equations of dynamics analytically, therefore numerical methods will be used. Development of these numerical methods is an important part of this work. An efficient simulation method that includes stochastic effects is the kinetic Monte Carlo algorithm, which is based on the chemical master equation.

Nonetheless this method includes a number of approximations and it is necessary to confirm its results by computer simulations at the microscopic level.

In the first part of Chap. II influence of fluctuations on a wave front (FKPP model) and a bistable system (Schlögl model) in which fast reactions cause perturbation of particle velocity distributions will be investigated. An important challenge is the derivation of a master equation that takes into account deviation of velocity distributions from the Maxwellian. Correctness of this perturbed master equation will be determined by comparison with microscopic simulations.

The second part of Chap. II is devoted to the characterization of the influence of internal fluctuations on formation of spatial structures, with emphasis on the development of shape during embryonic growth. A method for simulating formation of Turing patterns at the microscopic level will be presented. Then, this method will be used to explain results of experiments known from developmental biology to determine if the Turing mechanism is in agreement with them. Next, a simple reaction model will be developed for the "clock and wavefront" mechanism, which is one of the postulated mechanisms of formation of one-dimensional periodic structures in vertebrate embryos. Such model allows numerical simulations using the kinetic Monte Carlo algorithm and, in consequence, studies of the influence of fluctuations on the "clock and wavefront" mechanism.

Selected numerical methods for simulating reaction-diffusion systems

So far we have dealt with different levels of analytic description of chemical systems:

the deterministic, macroscopic one and the stochastic, mesoscopic one. In both cases, except for simple models, the equations for dynamics have no analytic solution. Therefore, to obtain information about evolution of the system, it is necessary to use numerical methods.

Computationally most efficient are algorithms for solving sets of macroscopic equations. They are presented in Sec. 4.1. These algorithms give us a deterministic trajectory of the system, which depends only on the initial conditions.

At the mesoscopic level, direct numerical integration of the master equation is highly inefficient due to a large number of degrees of freedom. Therefore, in order to solve the master equation, we use the kinetic Monte Carlo algorithm, described in Sec. 4.2.

Simulations using this algorithm give a single, stochastic trajectory of the system, which is equivalent to sampling from the probability distribution of system states. To get full information about this probability distribution it is necessary to repeat the simulation many times.

A separate group of numerical methods are microscopic simulations. They are not derived from analytic methods, but aim at simulating motion of all molecules in the system. These methods are the closest to real systems, but also the least efficient. Two microscopic methods, molecular dynamics and direct simulation Monte Carlo, will be presented in Sec. 4.3.

Numerical solution of macroscopic equations

There are many iterative methods for solving sets of first-order ordinary differential equations, such as macroscopic equations for a homogeneous system (I.9) [49]. Given initial conditions, iterations of the algorithm give state of the system at consecutive moments in time, which requires a time step ∆t to be assumed. Conceptually the simplest is the Euler algorithm. In it, the system evolves according to:

c i (t + ∆t) = c i (t) + f i (c,t)∆t. (I.63)
This method has serious drawbacks: it is numerically unstable and its global truncation error is high, proportional to ∆t. A more precise, and therefore allowing longer time steps, is the Runge-Kutta method [50]. In this method for a one-variable system, system state in consecutive time steps is calculated as:

c(t + ∆t) = c(t) + 1 6 ∆t(k 1 + 2k 2 + 2k 3 + k 4 ) (I.64) k 1 = f (c,t) (I.65) k 2 = f c + ∆t 2 k 1 ,t + ∆t 2 (I.66) k 3 = f c + ∆t 2 k 2 ,t + ∆t 2 (I.67) k 4 = f (c + ∆tk 3 ,t + ∆t) . (I.68)
The Runge-Kutta algorithm is numerically stable and its total truncation error is proportional to O (∆t 4 ).

In the case of reaction-diffusion systems, which are described by equations in the form (I.32), the above-mentioned methods can also be used, but first, these equations have to be transformed into ordinary differential equations [51]. One way such approximation can be done is by the finite difference method. It is based on discretization of space into cells of size ∆x. For a one-dimensional system the simplest approximation of a gradient between neighbouring cells is

∇c x + ∆x 2 = c(x + ∆x) -c(x) ∆x . (I.69)
Similarly, second derivative in a given cell can be calculated as

∆c (x) = c(x -∆x) + c(x + ∆x) -2c(x) ∆x 2 . (I.70)

Mesoscopic simulations -the kinetic Monte Carlo method

Based on the master equation (cf. Sec. 1.2) Gillespie developed an algorithm for simulating a stochastic trajectory of a system, called the kinetic Monte Carlo method [52,53].

This mesoscopic method has an advantage that, while requiring relatively low computational cost, it includes internal fluctuations without any additional assumptions about their probability distribution. On the other hand, it has all the limitations of master equation. It is applicable only to Markov processes and an expression for transition probabilities, which usually contain macroscopic parameters, has to be known.

In the case of inhomogeneous systems, state of a system is described by a vector of populations of species in each spatial cell. Cell size has to be carefully chosen. They have to be small enough to be homogeneous, but large enough for the transport between cells to be of diffusive nature. Moreover, simulation efficiency significantly decreases with shrinking of cells.

A simulation run is initiated by defining species populations in each cell. A single simulation step represents a single process taking place in the system. In reactiondiffusion systems, such process can be either a single chemical reaction or a diffusive jump of a molecule to a neighbouring cell. Algorithm step starts with calculating transition probabilities for all M processes ω i , i = 1,...,M and their sum ω T . Probability that no process will take place over time period τ is equal to P (τ ) = exp (-ω T τ ). So it provides a distribution of waiting time for any process to happen. Temporal length of a simulation step is sampled from this distribution and equal to

∆t = 1 ω T ln 1 r 2 , (I.71)
where r 2 is a random number with uniform distribution over (0,1]. The probability, that after time ∆t an i-th process will take place is equal to ω i /ω T . In order to select a process a number r 1 is sampled from a uniform distribution over (0,1]. In the given simulation step process number β will take place, for which

β-1 i=1 ω i < r 1 ω T ≤ β i=1 ω i . (I.72)
Last part of a simulation step is updating the state vector according to the selected process and increasing simulation time by ∆t.

Microscopic simulations

At a finer level a reaction-diffusion system is a mixture of particles that move around and collide with each other. Thus the most accurate simulation methods in the framework of classical mechanics rely on tracking positions and momenta of all particles, what we call microscopic simulations [54]. This approach is closer to reality and helps better understand phenomena that are difficult to study experimentally. Use of simpler mesoscopic or macroscopic methods can be justified by comparing them with microscopic simulations. The major drawback of microscopic algorithms is their large computational cost, which limits them to small systems and short simulation times.

Choice of a microscopic simulation method depends on the density of the system of interest. In dense media intermolecular collisions are frequent and correlations between positions and momenta of a collision pair are significant. Moreover, long-range intermolecular forces can be important. In that case a molecular dynamics (MD) method is used, in which Newton's equations of motion are used to update position and momenta of all the particles [55]. The MD method becomes inefficient for rarefied systems, in which mean free path is much larger than particle sizes. Then, the direct simulation Monte Carlo (DSMC) method is more adequate [18].

Another important issue is the method of simulating chemical reactions. The best way would be to include molecular orbitals of all particles and perform chemical reactions according to rules of quantum mechanics [56]. However, due to computational cost, such detailed treatment of reactions can be done only for systems of a few particles, while this work deals with collective behaviour of systems containing tens of thousands of particles.

From the global point of view, details of single reactions are not important, so simplified models of reactions are used. Moreover particles are considered to be spherical.

Direct simulation Monte Carlo method

The DSMC is a method of microscopic simulations of rarefied media proposed by Bird [18,57]. It can be seen as a method for numerically solving the nonlinear Boltzmann equation [58]. It was developed to model flows in high-atmosphere flights, but since then it has been used in other fields, like nanoscopic flows or chemical fluctuations. The method is based on the assumption that motion of molecules can be separated from their collisions. In each simulation step ∆t, particle positions are changed according to the equations of motion r(t + ∆t) = r(t) + v∆t. Given particle, while travelling the distance v∆t, may have collided with one of the surrounding particles. Hence the system is divided into spatial cells and we assume that only particles belonging to the same cell can collide. Introducing cells also accounts for spatial inhomogeneities in the system. The number of collisions, collision pairs and particle post-collision velocities are sampled from distributions known from the kinetic theory of gases. Separation of motion of particles from collisions can be used only if the time step ∆t is much shorter than the mean free time. Then, probability that a particle collides more than once in a single time step is negligible.

The average number of collisions per unit time within a cell is a product of the number of pairs of particles in the cell N (N -1)/2, where N is the number of particles in the cell, and the probability of collision of a pair per unit time. This probability of collision depends on the model of intermolecular interactions used. Here, we will focus on the hard spheres model, but there are more elaborate models fit for engineering applications [18,59]. A collision of hard spheres takes place when the distance between centres of particles is smaller than the sum of their radii, denoted by σ. In the frame of reference of one of the molecules of the collision pair, the other one combs through the volume V coll = πσ 2 v r ∆t per time step, where v r is the relative velocity between them. Since we assumed that cells are homogeneous, probability of collision of a given pair is equal to V coll /V , where , where R is a random number uniformly distributed over [0,1). This algorithm of selecting collision pairs is in agreement with binary collisions integral from the Boltzmann equation and is more efficient than explicitly considering probabilities of collisions between all pairs of molecules.

V is the cell volume. In each cell N (N -1)πσ 2 v max
Next, geometry of a collision has to be defined. All collision parameters are equally probable, because cells are homogeneous. Particles pre-collision velocities are denoted by v 1 i v 2 , while their relative velocity is v r = |v 1 -v 2 |. We select a frame of reference of the centre of mass of the pair, which travels with the velocity v m . In this frame, elastic scattering caused by collisions of hard spheres is isotropic, so probability distribution of the direction of the post-collision relative velocity vector v * r is spherically uniform. Coordinates of the post-collision relative velocity vector are [18] v * r = (cos(β)v r , sin(β) cos(ǫ)v r , sin(β) sin(ǫ)v r ) , (I. 73) where cos(β) and ǫ are sampled from a uniform distribution over [-1,1] and [0,2π), respectively. Particles post-collision velocities in the frame of reference of the system are

v * 1 = v m + v * r 2 (I.74) v * 2 = v m - v * r 2 . (I.75)
Collisional chemical reactions are treated as elastic collisions that result in a change of chemical identities of particles. A reaction takes place if energy of the collision is greater than the activation energy of the reaction. There may be additional, geometric constraints on the reaction, caused by internal structure of particles. To take this into account, the so-called steric factor is introduced. It is a value of probability in the range (0,1], that a given collision with enough energy will in fact result in a reaction. The first-order reactions are independent of collisions and their number per time step in each cell is sampled from the binomial distribution.

Molecular dynamics

The term molecular dynamics stands for computer simulations in which motion of all molecules of a mixture is simulated by changing their positions and momenta according to the Newton's laws of motion [55]. Therefore it is a deterministic method. A crucial matter is the choice of the intermolecular potential U (r). The total force not associated with chemical bonds or electrostatic interactions between ions is called the Van der Waals force. It includes the short-range repulsion due to Pauli's exclusion principle and longrange attraction between stable or induced dipoles. The most popular approximation of the Van der Waals force is the Lennard-Jones (LJ) potential [60]:

U LJ (r) = 4ǫ σ r 12 - σ r 6 , (I.76)
where r is the distance between centres of particles, ǫ is the depth of the potential minimum and σ is a distance, at which the potential is equal to zero. The long-range attraction is relatively weak and in some cases, e.g. at high temperature or in a system of low density, can be neglected. In that case the Weeks-Chandler-Andersen potential is used [61]:

U W CA (r) =      4ǫ ( σ r ) 12 -( σ r ) 6 + ǫ r < 2 1/6 σ 0 r ≥ 2 1/6 σ . (I.77)
The WCA potential is identical to the LJ potential, but cut at the position of its minimum and shifted by ǫ to remove discontinuity.

Force acting on the i-th particle due to interactions with the j-th particle is equal to

F ij = -d r U (|r|) r |r| , (I.78)
where r = r ir j . The net force acting on the j-th particle is F j = i =j F ji . In each step of the simulation net forces acting on all particles have to be calculated. It is part of the algorithm that requires greatest computational effort, proportional to the square of the total number of particles. To obtain linear scaling with system size, an approximation, in which the intermolecular potential is cut at a distance at which its gradient is negligible, is used. For the Lennard-Jones potential it is usually the distance r c = 2.5σ. Simulations that use the WCA potential are faster, because for this potential the cutoff is r c = 2 1/6 σ. After net forces are calculated, equations of motion are numerically integrated. Selected methods of numerical integration were presented in Sec. 4.1. In the simplest algorithm, particle momenta are changed by ∆v = F ∆t/m and their positions by ∆x = v∆t. Length of a time step ∆t is usually of a few femtoseconds.

If the goal is to simulate a canonical system, i.e. a system in thermal contact with its surroundings, it is necessary to include an algorithm for the thermostat. One of such algorithms is the Nosé-Hoover thermostat [62]. In it, phase space is extended by an additional variable η. Its equation of motion is

d t η = E k -gk B T, (I.79)
where E k is the total kinetic energy, g is the number of degrees of freedom and T is the temperature of the thermostat. Now, net force acting on a single particle contains an additional term representing exchange of energy with the thermostat. The net force acting on the j-th particle is

F j = - i =j d r U (|r|) r |r| - η Q p j , (I.80)
where Q is a parameter controlling strength of thermal coupling between the system and the surrounding, while p j is momentum of the j-th particle.

Effect of perturbation of particle velocity distribution on system dynamics

Rate of chemical reactions can be calculated from the frequency of collisions of particles, which is usually derived assuming a Maxwellian distribution of velocities of particles of each species [17,63]. In most cases this assumption is valid, but if the reaction is very fast, i.e. when frequency of reactive collisions is comparable to frequency of elastic collisions, distributions of particle velocities are significantly perturbed [64][65][66][67][68]. This effect is caused by the fact that particles with higher energy have higher chance of reacting than low-energetic ones. This has two reasons: one of them is that frequency of collisions increases with kinetic energy, the other is the effect of activation energy, which is the minimal energy of a collision pair required to react. In consequence, a shift towards higher velocities is seen in velocity distributions of products of the reaction. If the system is thermoneutral, i.e. the reaction has no thermal effect, total velocity distribution of all species is not perturbed and has the equilibrium form. Hence, there is a corresponding shift of distributions of substrates towards lower velocities.

The perturbation has an influence on reaction dynamics, which now is not described by the standard rate equations. The correct equations for dynamics contains some additional terms that are derived from the Boltzmann equation by keeping higher-order terms of the perturbative solution [69]. In general, there are corrections both to the rate of reaction and the diffusive flow.

Multistability

Perturbation of velocities distribution can affect number of stationary states of a system, what was shown on the Schlögl model of bistability (cf. Sec. 2.1.1). The perturbation causes a shift of the domain of bistability in the parameter space. It means that a system that, using traditional analysis, would be monostable, in reality is bistable.

Other types of transition were also observed: from bi-to monostability and between two monostable systems, but with significantly different positions of the stable state. Finally, it was proven that the perturbation may cause two stable states to move closer to each other, resulting in a decrease of the mean passage time by several orders of magnitude.

The perturbation increases with decreasing activation energy and increasing steric factor. Results obtained using a master equation that includes the perturbation, which has been derived in this work from the macroscopic equation [68], were compared with DSMC simulations. We found that the perturbed master equation is no longer valid for activation energies lower than 6k B T .

Wave front

Influence of perturbation of velocity distributions on the development of a wave front was investigated using the FKPP model. Since pushed wave fronts are sensitive to fluctuations, a mesoscopic approach was applied. Master equation that includes the perturbation of particle velocity distributions was derived from a perturbed macroscopic equation [70,71]. In the second-order approximation, perturbation of velocity distributions affects not only the rate of reactions, but also the diffusion of species [70].

We found that for fast reactions width and speed of the wave front can be up to 40% higher than the non-perturbed equations predict. This result was confirmed by microscopic simulations using the DSMC method.

Formation of one-dimensional spatial patterns

Due to coupling between chemical reactions and diffusion, which usually leads to higher homogeneity, stable spatial patterns may appear. This surprising ability of reaction-diffusion systems was first shown by Alan Turing in his work on morphogenesis, which is the development of shape in animal kingdom [43]. It is in the context of biology that this spatial structures are of key importance [26,72,73]. Nevertheless, they were also observed in chemical reactors, what directly proved their existence [74,75].

Simulating formation of Turing patterns with the DSMC method

The major difficulty in simulating formation of Turing patterns at the microscopic level lies in obtaining significantly different diffusion coefficients for two species. In a mixture of hard spheres, particle mobilities depend on their masses and sizes. In order to get two species diffusing with different rates, the mixture has to contain at least three types of particles. Two of them react with each other, while the third one acts as solvent.

Contrary to the analysis of Turing patterns done in Sec. 2.2.2, diffusivities in such system are not constant, but depend on the local composition of the mixture [17]. However, they can be considered constant in the area of phase space around the state, from which the pattern develops, if population of solvent particles is much larger than that of other species. Then, majority of collisions of reactive species are with solvent molecules, while collision between reactive species have a negligible impact on their diffusivities. In the simulation algorithm, solvent particles play an additional role of the reservoir, from which reactive particles are created at a constant rate, what keeps the system out of equilibrium. Computational efficiency is greatly increased if collisions between solvent particles, that have no effect on the structure in concentrations of other species, are not simulated. Using the developed method we observed that internal fluctuations may cause transitions between structures of different wavelengths.

Somitogenesis

An animal embryo grows from a single fertilised cell, from which further, similar cells form that can differentiate into more specialised cells. Development of an embryo is a process of coordinated differentiation of groups of cells, which creates the complex structure of an animal body. How do cells obtain information about their positions is not fully understood [76]. It is known that some substances, called morphogens, affect differentiation path of cells. Therefore spontaneous formation of spatial order in reaction-diffusion systems, e.g. Turing patterns could also explain many characteristic features of somitogenesis [81][82][83].

Somites are formed one after another starting from the head and going towards the tail. They are laid down symmetrically on both sides of the neural tube. Simultaneously, the embryo elongates through formation of new cells at the tail bud. The layer of undifferentiated cells between somites and the tail is called the presomitic mesoderm (PSM).

The moving boundary between cells with already defined differentiation path (although they may have not yet differentiated into somites) and the presomitic mesoderm is called the determination front.

The exact, molecular mechanism of somitogenesis is not known and its details can be specific to a given species. However, it is believed that a general mechanism of formation of somites is common to all vertebrates and it is possible to find key morphogens involved in it. Propagation of the determination front is attributed to the retinoic acid (RA), which is found in the area already divided into somites, and the fibroblast growth factor (FGF), which has high concentration in the presomitic mesoderm [84,85]. Another important observation is that in embryos of many vertebrates oscillations of expression of some genes are present in the PSM [86,87].

A theoretical analysis of perturbation of somitogenesis

In order to experimentally determine the mechanism of somitogenesis, formation of somites is artificially perturbed. One way such perturbation can be imposed is by grafting a source of morphogens close to the anterior-posterior axis during somitogenesis [88].

In result, somite close to the source is enlarged, while those anterior of it appear smaller. This result was so far considered to support the "clock and wavefront" mechanism (cf. Sec. 6.2.2). However in our work, using the developed method of simulating Turing patterns, we show that the Turing mechanism can also account for the experimental observations. We analyse in detail the influence of position, size and moment of introduction of the bead on the resulting phenotype. Moreover, we study a case in which there is a local shortage of morphogen in the embryo, which opens new experimental opportunities of identifying the mechanism of somitogenesis.

Reaction-diffusion scheme for the "clock and wavefront" model

The first model of the "clock and wavefront" type was proposed by Cooke and Zeeman to explain somitogenesis [78]. Many similar models were developed since then. Their common feature is that the structure of somites is formed due to an interaction between an intermolecular oscillator (clock) and a propagating wave front. In the original model it is assumed, that synchronised oscillators work in the undifferentiated cells in the posterior part of the embryo and that the wave front propagating from the anterior towards the posterior end stops the oscillators in their current phase. Consequently a periodic structure forms, wavelength of which is equal to the distance travelled by the wave front during one period of the intermolecular oscillations. This model gained a lot of attention after gene expression oscillations were observed in the posterior parts of embryos of some species. So far, mathematical representation of the model was done using macroscopic equations that have no microscopic basis. In our paper we present for the first time a scheme of reactions that, together with diffusion, results in the behaviour proposed by Cooke and Zeeman. The criterion of choice of the scheme was its simplicity -we were looking for a minimal set of reactions that would give the required behaviour. The scheme is built from two connected parts -a model of an antagonistic wave front, derived from the Schlögl model, and a model of chemical oscillator, inspired by the Brusselator. Species involved in the wave front are assumed to catalyse some of the reactions of the oscillator. Thanks to this coupling, there is a qualitative change in the dynamics of the oscillator at the moment of front passage.

Ahead of the front, the oscillator model has a stable limit cycle attractor, while after its passage it has two stable stationary states. Moreover, the separatrix between basins of attraction of these states cuts through the limit cycle. In consequence, after front passage, cells are in one of two stable states, depending on their phase of oscillations at the moment of front arrival.

Based on the developed reaction-diffusion model, simulations using the kinetic Monte Carlo algorithm were performed. They allowed us to rigorously study the effect of internal fluctuations on the formation of the structure. These simulations showed that to form a regular structure it is necessary that intermolecular oscillators are coupled with each other. If there is no coupling between them, oscillators gradually go out of phase and the final structure behind wave front is irregular. An extended model was proposed, in which neighbouring oscillators are coupled through diffusion.

III. Conclusions

An in-depth study of stochastic effects in some generic models of nonlinear phenomena was conducted. The analysis was done on three levels of accuracy, each of them revealing new effects. It can be concluded that the mesoscopic description based on the master equation is sufficiently accurate, while remaining computationally efficient. A method of extending this description to systems in which fast reactions perturb the particle velocity distributions (an effect best studied at the microscopic level) was presented.

The range of validity of the mesoscopic description was determined by comparison with microscopic simulations, which are associated with longer computation times and may be less appropriate in some fields, such as biology. On the other hand, the deterministic, macroscopic approach proved to be oversimplified to account for the effects discussed, since the phenomena of interest are greatly affected by fluctuations.

New techniques of simulations at the microscopic level were developed to model the formation of Turing patterns, which are of particular importance in biology. My studies show that the Turing mechanism should remain a possible explanation of axis segmentation in vertebrae embryos, alongside the "clock and wavefront" model. Ideas for further experiments that would help determine the correct mechanism were presented. Also, I successfully formulated the "clock and wavefront" model in the framework of reactiondiffusion systems, which offers a microscopic picture of the phenomenon, that could be useful to engineer an artificial system with analogous complex behaviour. The study also confirmed the importance of coupling of oscillators in this model. It is worth to note that the minimal reactions scheme for the "clock and wavefront" model is much more complex than the one necessary for the formation of Turing patterns.

Instytut Chemii Fizycznej

Polskiej Akademii Nauk 

Ujęcie makroskopowe

W ujęciu makroskopowym dynamikę układu składającego się z N reagentów opisuje układ N równań różniczkowych typu [12]:

∂ t c i (x,t) = f i (c 1 ,c 2 ,...,c N ) -∇ • j i , (I.1)
gdzie c i (x,t) to stężenie i-tego reagenta w punkcie x w czasie t, j i to strumień dyfuzji, a f i (c 1 ,c 2 ,...,c N ) = j a ij g j (c 1 ,c 2 ,...,c N ) to sumaryczne tempo zmiany stężenia i-tego składnika na skutek reakcji chemicznych [13,14]. Współczynnik a ij równy jest zmianie ilości cząsteczek i-tego składnika spowodowanej przez pojedynczą j-tą reakcję. W przypadku reakcji elementarnej, to jest zachodzącej bez produktów pośrednich, szybkość reakcji g j można policzyć z równania kinetycznego reakcji chemicznej [15,16]:

g j = k j c s 1j 1 c s 2j 2 ...c s N j N , (I.2)
gdzie k j oznacza stałą szybkości reakcji, a s ij to współczynnik stechiometryczny przy i-tym składniku w mechanizmie reakcji j. Przykładowo, dla reakcji A + B

k 1 --⇀ ↽ -- k -1 2A, otrzymamy f A = -f B = k 1 c A c B -k -1 c 2
A . W ogólności, strumień dyfuzji i-tego składnika j i zależy od niejednorodności przestrzennych w stężeniach wszystkich składników mieszaniny [17]. Jednak w wielu przy-padkach dyfuzja krzyżowa jest zaniedbywalna i można przyjąć założenie, że dyfuzja danego składnika jest proporcjonalna do gradientu jego stężenia. Wyraża się to prawem Ficka j i = -D i ∇c i , gdzie D i to współczynnik dyfuzji i-tego składnika. Liniowa zależność strumienia dyfuzji od gradientu stężenia jest poprawna, gdy liczba Knudsena, wyrażająca stosunek średniej drogi swobodnej cząsteczek do długości charakterystycznej niejednorodności stężeń, jest mniejsza od 0.1 [18]. Jeżeli współczynnik dyfuzji jest stały, to równanie (I.1) przyjmuje postać:

∂ t c i (x,t) = f i (c 1 ,c 2 ,...,c N ) + D i ∆c i . (I.3)
Poza najprostszymi przypadkami, układów równań typu (I.3) nie da się rozwiązać analitycznie i do ich rozwiązywania stosuje się algorytmy numeryczne opisane w rozdziale 4.1.

Ujęcie mezoskopowe -równanie master

Naturalnymi zmiennymi opisującymi stan układu chemicznego są populacje (czyli liczba cząsteczek) poszczególnych składników, które oznaczamy jako N i (t). Na razie ograniczymy się do układów jednorodnych. Stan układu opisuje wektor N , którego składowe to populacje wszystkich składników. Jako że na tym poziomie opisu chcemy uwzględnić fluktuacje wewnętrzne układu, interesuje nas ewolucja w czasie rozkładu prawdopodobieństwa stanów P (N ,t). Reakcje powodują skokowe zmiany w populacjach, więc dynamikę układu chemicznego można opisywać poprzez przypływ i wypływ (birth and death) prawdopodobieństwa. Zakładając, że prawdopodobieństwa następujących po sobie reakcji są od siebie niezależne, tj. są procesami Markowa, można wyprowadzić równanie na pochodną d t P (N ,t), zwane równaniem master [5,19,20]:

d t P (N ,t) = j
w j (Nr j )P (Nr j ,t)j w j (N )P (N ,t). (I.4)

Wektor r j oznacza zmianę wektora stanu pod wpływem j-tej reakcji. Jako w j (N ) oznaczamy prawdopodobieństwo przejścia ze stanu N do stanu N + r j na skutek j-tej reakcji. Pierwsza suma w tym równaniu odpowiada przypływowi prawdopodobieństwa do stanu N , natomiast druga suma to wypływ prawdopodobieństwa z tego stanu do innych stanów. Sumy te obejmują wszystkie reakcje.

Dla reakcji drugiego rzędu prawdopodobieństwa przejść są proporcjonalne do liczby par cząsteczek, które mogą ulec kolizji [5]. Na przykład dla modelu A + B

k 1 --⇀ ↽ -- k -1
2A, przytoczonego w rozdziale. 1.1, wynoszą one:

w 1 (N A ,N B ) = k m 1 N A N B (I.5) w -1 (N A ) = k m -1 2 N A (N A -1). (I.6)
Mają one postać podobną do równania na szybkość reakcji, jednak stałe

k m 1 i k m -1
różnią się od makroskopowych stałych szybkości reakcji, ponieważ tutaj nie odnoszą się do stężeń, lecz do populacji. Ponadto w wyrażeniu na w -1 pojawia się czynnik

kombinatoryczny 1 2 . Mamy więc k m 1 = k 1 /V i k m -1 = 2k -1 /V
, gdzie V oznacza objętość układu. Równanie master jest równaniem bardziej fundamentalnym od równania reakcji-dyfuzji, które można z niego wyprowadzić [19]. Deterministyczna, makroskopowa zmienna określającą liczbę cząsteczek i-tego rodzaju definiowana jest jako średnia wartość N i : N i = N N i P (N ). Jeśli użyjemy postaci całkowej równania master, jako pierwsze przybliżenie d t N i otrzymamy równanie makroskopowe. Pierwszy człon tego równania oznacza przypływ i wypływ prawdopodobieństwa na skutek reakcji, natomiast drugi -na skutek dyfuzji. Prawdopodobieństwa przejścia dla dyfuzji zależą od wybranego podziału przestrzeni na komórki. W ogólności można je liczyć na podstawie prawa Ficka stosując metodę różnic skończonych [21]. W przypadku jednowymiarowym, gdy wszystkie komórki mają takie same rozmiary, dywergencję z gradientu populacji w k-tej komórce można liczyć jako

d t P ({N }) = k j [w j (N k -r j )P (N k -r j ) -w j (N k )P ] + k l [d k→l (N k + 1)P (N k + 1, N l -1) -d k→l (N k )P ] , ( 
(N k-1 + N k+1 -2N k )/(∆x) 2 ,
gdzie ∆x to wielkość komórki [5]. Odpowiadające temu wyrażeniu prawdopodobieństwo przejścia dyfuzyjnego wynosi

d k→l (N k ) = DN k ∆x 2 δ k±1,l , (I.8)
gdzie δ i,j oznacza deltę Kroneckera.

Równanie master jest podejściem pośrednim pomiędzy makroskopowymi równaniami reakcji-dyfuzji a opisem mikroskopowym. Poza niektórymi prostymi przypadkami, nie ma ono analitycznego rozwiązania. Jednakże na jego podstawie możliwe są symulacje dynamiki układu za pomocą metody kinetic Monte Carlo (KMC), opisanej w rozdziale. 4.2.

Nieliniowe układy dynamiczne

Rozdział ten zawiera wprowadzenie do nieliniowych układów dynamicznych [22][23][24].

Różnorodność 

Układy jednorodne

d t c 1 = f 1 (c)
. . . (I.9)

d t c N = f N (c).
Stany stacjonarne można znaleźć rozwiązując układ równań d t c = 0. Stabilność stanu stacjonarnego c s określamy poprzez rozwinięcie równań z układu (I.9) w szereg Taylora wokół tego stanu [23,24]. Ponieważ interesują nas niewielkie wychylenia u = cc s , zachowujemy jedynie człony liniowe. Uwzględniając f i (c s ) = 0, otrzymujemy: [26,27]:

d t u T = J s u T , J =       ∂f 1 ∂c 1 . . .
• λ 1 ∈ R, λ 2 ∈ R, λ 1 = λ 2 , λ 1 < 0, λ 2 < 0
Stabilny węzeł -układ w sposób eksponencjalny w czasie zmierza do stanu stacjonarnego wzdłuż wektorów własnych Jakobianu. Rysunek I.1(a) przedstawia trajektorie układu dla przypadku λ 2 < λ 1 , w którym w pobliżu stanu stacjonarnego układ porusza się wzdłuż wektora własnego v 1 .

•

λ 1 ∈ R, λ 2 ∈ R, λ 1 = λ 2 , λ 1 > 0, λ 2 > 0
Niestabilny węzeł -układ w sposób eksponencjalny w czasie opuszcza stan stacjonarny wzdłuż wektorów własnych Jakobianu.

• λ 1 ∈ R, λ 2 ∈ R, λ 1 = λ 2 , λ 1 < 0, λ 2 > 0
Punkt siodłowy -niestabilny stanu stacjonarny, pokazany na Rys. I.1(b). Jedynie niewielkie wychylenia wzdłuż wektora v 1 będą zanikały. Wzdłuż wektorów własnych przez punkt siodłowy przechodzą separatrysy.

• λ 1 ∈ C, λ 2 ∈ C, λ 1 , λ 2 = α ± iβ, α < 0
Stabilne ognisko -rozwiązanie układu równań (I.10) zawiera człony proporcjonalne do exp[±iβt], więc układ zmierza do stanu stacjonarnego w sposób oscylacyjny, co jest pokazane na Rys. I.1(c).

• λ 1 ∈ C, λ 2 ∈ C, λ 1 , λ 2 = α ± iβ, α > 0
Niestabilne ognisko -układ w sposób oscylacyjny opuszcza stan stacjonarny. Prostym modelem układu bistabilnego, to jest posiadającego dwa stabilne stany stacjonarne, jest model Schlögla [28,29]:

• λ 1 ∈ C, λ 2 ∈ C, λ 1 , λ 2 = ±iβ
2X + A k 1 --⇀ ↽ -- k -1 3X (I.11) X k 2 --⇀ ↽ -- k -2
B, (I.12) gdzie stężenia składników A i B przyjmuje się za stałe. Założenie takie można przyjąć gdy te stężenia są znacznie większe od stężenia składnika X, lub gdy są one utrzymywane na stałym poziomie poprzez wymianę z zewnętrznymi rezerwuarami, jak w przypadku reaktora typu CSTR (continuous stirred tank reactor). Dynamikę układu opisuje jedno równanie, mające postać wielomianu trzeciego stopnia: Stacjonarny rozkład prawdopodobieństwa można policzyć z warunku na zanik wypadkowych strumieni prawdopodobieństwa pomiędzy stanami. Załóżmy, że możliwe są przejścia tylko między stanami, których populacje różnią się o 1, tak jak w modelu (I.11,I.12).

d t X = -k -1 X 3 + k 1 AX 2 -k 2 X + k -2
Przyjmując dowolną wartość P (0), prawdopodobieństwa innych stanów wyrażają się przez wzór iteracyjny w + (N X )P (N X ) = w -(N X + 1)P (N X + 1), gdzie w + to prawdopodobieństwo przejścia do stanu o populacji większej o jeden, a w -to prawdopodobieństwo przejścia w odwrotną stronę.

Układ może pozostawać przez długi czas w okolicy któregoś z maksimów, jednak odpowiednio duże fluktuację mogą spowodować przejście w okolice drugiego stanu stabilnego.

Średni czas takiego przejścia można obliczyć wykorzystując prawdopodobieństwa przejścia z równania master (I.4) [19]: Jednym ze sposobów pojawienia się stabilnego cyklu granicznego jest nadkrytyczna bifurkacja Hopfa [31]. Polega ona na tym, że przy zmianie wartości parametrów następuje przejście od wartości własnej Jakobianu o ujemnej części rzeczywistej, przez zero, do Reλ > 0, przy czym części urojone wartości własnych są niezerowe, czyli jest to przejście od stabilnego do niestabilnego ogniska.

τ 0→2 = N X2 -1 ν=N X0 ν µ=0 1 w + (µ) ν ξ=µ+1 w -(ξ) w + (ξ) (I.14) τ 2→0 = N X2 -1 ν=N X0 µ≥ν 1 w -(µ) µ ξ=ν+1 w + (ξ) w -(
Rozważamy układ dwóch zmiennych:

d t c 1 = f (c 1 ,c 2 ), d t c 2 = g(c 1 ,c 2 ). (I.16)
Przyjmijmy, że układ ten posiada tylko jeden stan stacjonarny (c s 1 ,c s 2 ). Dla ułatwienia analizy jego stabilności wprowadzamy nowe zmienne x = c 1 -c s 1 , y = c 2 -c s 2 , dzięki czemu badany stan stacjonarny znajduje się w punkcie x = 0,y = 0. Linearyzując równania (I.16) wokół tego stanu otrzymujemy:

d t x = f x | (0,0) x + f y | (0,0) y (I.17) d t y = g x | (0,0) x + g y | (0,0) y. (I.18)
Stabilność tego stanu określają wartości własne Jakobianu powyższego układu. Można je policzyć na podstawie wielomianu charakterystycznego Jakobianu w punkcie stacjonarnym:

det (J s -λI) = λ 2 -λtr(J s ) + det(J s ) = 0. (I.19)
Wartości własne spełniające to równanie wynoszą:

λ 1,2 = tr(J s ) ± tr 2 (J s ) -4 det(J s ) 2 . (I.20)
Stan (0,0) będzie więc niestabilnym ogniskiem pod warunkiem, że: tr(J s ) > 0, tr 2 (J s ) < 4 det(J s ). A

k 1 -→ X (I.22) B + X k 2 -→ Y + D (I.23) 2X + Y k 3 -→ 3X (I.24) X k 4 -→ E. (I.25)
Jest to model reakcji w układzie otwartym, w którym następuje produkcja składnika E ze składnika A oraz składnika D ze składnika B. Dla wygody stosujemy następującą zamianę zmiennych

X ′ = k 3 k 4 X, Y ′ = k 3 k 4 Y, t ′ = k 4 t (I.26) A ′ = k 2 1 k 3 k 3 4 A, B ′ = k 2 k 4 B. (I.27)
W nowych zmiennych, pomijając primy w oznaczeniach, równania kinetyczne mają postać

d t X = X 2 Y -(B + 1)X + A (I.28) d t Y = -X 2 Y + BX. (I.29)
Model ten posiada tylko jeden stan stacjonarny X 0 = A, Y 0 = B A . Jakobian w tym punkcie ma postać:

J s =    B -1 A 2 -B -A 2    . (I.30)
Na podstawie warunków (I.21) wiemy, że nadkrytyczna bifurkacja Hopfa ma miejsce w punkcie:

B = A 2 + 1. (I.31)
Cykl graniczny jest to atraktor strukturalnie stabilny. Oznacza to, że wychylenia od jego trajektorii będą zanikały, co go odróżnia od oscylacji w układzie zachowawczym, w którym trajektoria zależy od warunków początkowych [26,33,34]. Tak więc cykl gra- 

Układy reakcja-dyfuzja

W wielu rzeczywistych układach mogą wystąpić przestrzenne różnice w stężeniach reagentów. W takim przypadku w równaniach dynamicznych uwzględniany jest dyfuzyjny transport składników. Dynamikę opisuje układ równań typu:

∂ t c i (x,t) = f i (c 1 ,c 2 ,...,c N ) + D i ∆c i , (I.32)
gdzie D i to współczynnik dyfuzji i-tego składnika.

Fronty falowe

Rozważmy niejednorodny układ chemiczny, którego dynamika zezwala na przebywanie przez długi okres czasu w jednym z dwóch stanów stacjonarnych. Załóżmy, że część układu znajduje się w jednym, a reszta w drugim z nich. Wówczas na granicy między tymi obszarami nastąpi transport dyfuzyjny składników, który, w połączeniu z reakcjami chemicznymi, spowoduje powiększanie się jednego z tych obszarów kosztem drugiego. Rozchodzący się gradient stężeń nazywany jest frontem falowym. W najprostszym przypadku jednowymiarowego układu z jedną zmienną dynamikę opisuje równanie: 

∂ t c(x,t) = f (c) + D∂ xx c. (I.
lim z→∞ c(z) = 0, lim z→-∞ c(z) = 1. (I.36)
Fronty rozchodzące się pomiędzy stanem quasi-stabilnym a stabilnym dzielimy na dwa rodzaje: fronty ciągnięte (pulled fronts) i fronty pchane (pushed fronts) [36]. Fronty ciągnięte charakteryzują się tym, że pochodna f ′ (c) ma najwyższą wartość w stanie stacjonarnym c = 0 znajdującym się przed czołem frontu. Wówczas prędkość propagacji frontu zdeterminowana jest poprzez dynamikę reakcji na jego czele, gdzie f (c) ≈ 0, a więc można zastosować przybliżenie liniowe funkcji f (c), co daje:

Dd zz c + vd z c + f ′ (0)c = 0. (I.37)
To równanie posiada realistyczne rozwiązanie w postaci funkcji nieujemnej pod warunkiem, że:

v ≥ v c = 2 Df ′ (0). (I.38)
Tak więc prędkość stabilnych frontów ciągniętych jest ograniczona od dołu przez wartość krytyczną v c . Szerokość frontu definiuje się jako wartość bezwzględną odwrotności jego nachylenia w punkcie przegięcia. Przyjmując, że punkt przegięcia spełnia c = 1/2, szerokość wyraża się jako w = v/f (0.5), jest więc proporcjonalna do jego prędkości.

Ponieważ szybkość nie jest jednoznacznie określona, zależy ona od warunków początkowych: dla frontów węższych od v c /f (0.5) prędkość będzie dążyła do v c , natomiast fronty szersze będą propagowały odpowiednio szybciej [37][38][39].

Ważnym przykładem frontu ciągniętego jest model Fisher-Kolmogorov-Petrowsky-Piskunow (FKPP) [40,41]. Jest to ogólny model mający zastosowanie do badania wielu zjawisk, takich jak rozprzestrzenianie się korzystnego genu w populacji, czy rozchodzenie się chorób zakaźnych. W kontekście chemicznym oznacza on rozchodzenie się frontu na

skutek reakcji A + B k - → 2A przy zachowaniu całkowitej liczby cząsteczek c A + c B = 1.
Układ opisywany jest jedną zmienną c = c A o dynamice: 

f (c) = kc(1 -c). (I.
v p > v c .
Wiele modeli reakcja-dyfuzja opisujących rzeczywiste układy posiada więcej niż jeden stabilny stan stacjonarny i często pojawia się zagadnienie rozchodzenia się frontu pomiędzy dwoma takimi stanami. Ten przypadek jest jakościowo różny od dotychczas rozważanych, ponieważ pomiędzy dwoma stanami stabilnymi zawsze znajduje się jeden stan niestabilny, który rozdziela ich obszary przyciągania [42]. Najprostszy model takiego frontu, zwanego trigger front, opisuje równanie f (c) posiadające trzy miejsca zerowe c 0 < c 1 < c 2 . Szybkość jego propagacji jest ściśle określona i wynosi 

v = - c 2 c 0 f (c)dc
v = k -1 D/2(c 0 + c 2 -2c 1 )
, gdzie D to współczynnik dyfuzji [26]. Kierunek propagacji takiego frontu jest zdefiniowany przez względne prawdopodobieństwo stanów stacjonarnych. Przykładowo, dla rozkładu z Rys. I.3, front będzie przemieszczał się od bardziej prawdopodobnego stanu N X2 w kierunku stanu N X0 .

Struktury Turinga

Szukając możliwych mechanizmów morfogenezy, t.j. powstawania kształtów w królestwie zwierząt, Turing pokazał, że w odpowiednich warunkach reakcje chemiczne w połączeniu z dyfuzją mogą powodować spontaniczne powstawanie stabilnych struktur przestrzennych w stężeniach reagentów [26,43]. Struktury te mogą powstawać w układach o co najmniej dwóch składnikach i kluczowa jest różnica w mobilnościach tych dwóch rodzajów cząsteczek. Szerszą teorię tzw. struktur dyssypatywnych, t.j. powstających w układach podtrzymywanych z dala od równowagi poprzez przepływy termodynamiczne, rozwinął Prigogine [5,44,45].

Powstawanie struktur Turinga przedstawimy na przykładzie jednowymiarowego układu opisywanego dwiema zmiennymi u i v.

∂ t u(x,t) = f (u,v) + D u ∂ xx u (I.41) ∂ t v(x,t) = g(u,v) + D v ∂ xx v. (I.42)
Zakładamy, że punkt (0,0) to stan stacjonarny układu, czyli f (0,0) = g(0,0) = 0. Skalę odległości dobieramy tak, żeby uzyskać D u = 1. Stosunek stałych dyfuzji oznaczamy jako 

d = D v /D u .
trJ | (0,0) = (f u + g v )| (0,0) < 0, |J | = (f u g v -f v g u )| (0,0) > 0. (I.43)
Następnie przeprowadzimy analizę stabilności pełnego układu, uwzględniając dyfuzję.

Przedstawmy rozwiązania w postaci kombinacji liniowej członów harmonicznych

u(x,t) = i=∞ i=1 a i exp(λ i t + ik i x) (I.44) v(x,t) = i=∞ i=1 b i exp(λ i t + ik i x). (I.45)
Po linearyzacji funkcji f i g w okolicy punktu (0,0) i z uwzględnieniem (I.44,I.45) układ równań (I.41,I.42) przybiera postać

λ i u = f u (0,0)u + f v (0,0)v -k 2 i u (I.46) λ i v = g u (0,0)u + g v (0,0)v -k 2 i dv (I.47)
dla każdej liczby falowej k i . Współczynniki λ i są to wartości własne macierzy

   f u -k 2 i f v g u g v -k 2 i d    (0,0)
, (I.48) czyli spełniają:

λ 2 i + λ i k 2 i (1 + d) -(f u (0,0) + g v (0,0)) + h(k 2 i ) = 0, (I.49) h(k 2 i ) = dk 4 i -(f u (0,0)d + g v (0,0))k 2 i + f u (0,0)g v (0,0) -f v (0,0)g u (0,0). (I.50)
Każdej liczbie falowej odpowiadają dwie wartości własne λ ± i . Jeśli część rzeczywista którejś z nich jest dodatnia, stan stacjonarny jest niestabilny wobec perturbacji przestrzennych o danej liczbie falowej i w układzie powstanie stabilna struktura przestrzenna.

Może mieć to miejsce, albo gdy współczynnik przy pierwszej potędze λ i w powyższym równaniu jest ujemny, albo gdy h(k 2 i ) < 0. Pierwsza z tych możliwości nie może być spełniona, ponieważ z warunków (I.43) wynika f u (0,0) + g v (0,0) < 0.

(I.51) Druga z tych możliwości, czyli h(k 2 i ) < 0 dla jakiegoś k i > 0, może być spełniona tylko, gdy współczynnik stojący przy k 2 i jest ujemny, ponieważ f u (0,0)g v (0,0)-f v (0,0)g u (0,0) > 0 zgodnie z drugim warunkiem (I.43), czyli: warunkiem na powstanie struktury Turinga jest więc f u (0,0) > 0, g v (0,0) < 0. Z warunku

f u (0,0)d + g v (0,
(f u g v -f v g u )| (0,0) > 0 (patrz równanie (I. 43 
)) dodatkowo wynika, że f v (0,0)g u (0,0) < 0, co pozostawia dwie możliwości wyboru znaku tych pochodnych, które dają jakościowo różne struktury. Jeśli f v (0,0) > 0, g u (0,0) < 0, harmoniczne struktury obydwu składników będą w przeciwfazie, czyli minimum stężenia jednego z nich będzie odpowiadało maksimum drugiego. Natomiast dla f v (0,0) < 0, g u (0,0) > 0 struktury te będą w fazie [26].

Warunek (I.52) jest konieczny, ale niewystarczający dla uzyskania h(k 2 i ) < 0. Warunkiem wystarczającym jest h(k 2 m ) < 0, gdzie k 2 m to położenie minimum funkcji h(k 2 i ):

k 2 m = 1 2 f u + g v d , (I.53)
z czego wynika:

(df u (0,0) + g v (0,0)) 2 4d > f u (0,0)g v (0,0) -f v (0,0)g u (0,0). (I.54)

Z powyższej nierówności można wyprowadzić warunek na minimalny stosunek współczynników dyfuzji dla określonych wartości stałych reakcji. Ten krytyczny stosunek oznaczymy d c i jest on dany jako rozwiązanie równania Rysunek I.7.: Przykładowa zależność dyspersyjna dla modelu Graya-Scotta (linia przerywana). Pionowe linie kropkowane odpowiadają modom dozwolonym przez warunki brzegowe.

d 2 c f 2 u (0,0) + 2(2f v (0,0)g u (0,0) -f u (0,0)g v (0,0))d c + g 2 v (0,0) = 0. (I.
J s =    φ + κ X 2 + -2(φ + κ) -X 2 + -φ    . (I.61)
Istotne znaczenie ma dobór warunków brzegowych, tutaj przyjmijmy warunki Neumanna [48] ∂ 

x u(0,t) = ∂ x v(0,t) = ∂ x u(L,t) = ∂ x v(L,t) =

Cel pracy

Efekty charakterystyczne dla nieliniowych układów dynamicznych z dala od równowagi, które zostały przedstawione w rozdziale. 

U W CA (r) =      4ǫ ( σ r ) 12 -( σ r ) 6 + ǫ r < 2 1/6 σ 0 r ≥ 2
F j = - i =j d r U (|r|) r |r| - η Q p j , (I.79)
gdzie Q to parametr regulujący siłę sprzężenia między termostatem a układem, a p j to pęd j-tej cząsteczki. 

Front falowy

Wpływ zaburzenia rozkładu prędkości na rozwój frontu falowego zbadano na modelu FKPP. Do opisu zastosowano równanie master, jako że fronty pchane są czułe na fluktuacje. Postać równania master zawierającego zaburzenie rozkładu prędkości wyprowadzono na podstawie makroskopowego równania uwzględniającego to zaburzenie [70,71].

W drugim przybliżeniu zaburzenie rozkładu prędkości wpływa nie tylko na tempo reakcji chemicznych, ale też na dyfuzję składników [70].

Dla szybkich reakcji, szybkość i szerokość frontu może być nawet do 40% większa od tej przewidzianej nie uwzględniając zaburzenia rozkładu prędkości. Rezultat ten potwierdzono za pomocą symulacji mikroskopowych metodą DSMC.

Powstawanie struktur przestrzennych w układach jednowymiarowych

Dzięki odpowiedniemu sprzężeniu reakcji chemicznych z dyfuzją, która zazwyczaj prowadzi do ujednolicania układu, powstawać mogą stabilne struktury przestrzenne. Ta zaskakująca zdolność układów reakcja-dyfuzja została pierwszy raz wykazana przez Turinga w jego przełomowej pracy dotyczącej morfogenezy [43], czyli powstawania kształtu w królestwie zwierząt. To właśnie w kontekście biologii rozwoju te struktury przestrzenne mają kluczowe znaczenie [26,72,73], jednak są one też obserwowane w reaktorach chemicznych, dzięki czemu bezpośrednio dowiedziono ich istnienia [74,75]. Somity powstają jeden za drugim od strony głowy w kierunku ogona zarodka symetrycznie po obu stronach cewy nerwowej. Jednocześnie embrion ulega wydłużeniu poprzez powstawanie nowych komórek na końcu jego ogona. Warstwa jeszcze niezróżnicowanych komórek pomiędzy somitami a ogonem zarodka nazywana jest mezodermą przyosiową (paraxial mesoderm). Przesuwająca się granica pomiędzy komórkami, których ścieżka różnicowania jest już określona (chociaż nie wytworzyła się jeszcze struktura somitów), a mezodermą przyosiową, nazywa się frontem różnicowania (determination front). Dokładny, molekularny mechanizm somitogenezy nie jest znany, a jego szczegóły mogą być różne w zależności od gatunku. Pomimo tego uważa się, że ogólny mechanizm powstawania somitów jest wspólny dla wszystkich kręgowców i możliwe jest wyodrębnienie kluczowych dla niego morfogenów. Przemieszczanie się frontu różnicowania związane jest z kwasem retinowym, którego wysokie stężenie obserwuje się w obszarze zróżnicowanym na somity, i z czynnikiem wzrostu fibroblastów (FGF), który znajduje się w mezodermie przyosiowej [84,85]. Ponadto w zarodkach wielu kręgowców obserwowane są oscylacje ekspresji niektórych genów w komórkach mezodermy przyosiowej [86,87].

Teoretyczna analiza zaburzenia rozwoju kręgosłupa

W celu ustalenia mechanizmu somitogenezy przeprowadza się eksperymenty polegające na zaburzaniu rozwoju kręgosłupa w zarodku. Jednym z takich zaburzeń jest wszczepienie dodatkowego źródła morfogenu w okolice struny grzbietowej w trakcie powstawania somitów [88]. Obserwuje się wówczas powstanie powiększonego somitu na wysokości zaburzenia, natomiast somity pomiędzy źródłem a głową zarodka zostają pomniejszone. Wynik ten dotychczas interpretowano jako argument za poprawnością mechanizmu "clock and wavefront" (patrz rozdział 6.2.2). W naszej pracy, za pomocą wcześniej opracowanej techniki mikroskopowych symulacji struktur Turinga, pokazujemy, że mechanizm Turinga też może wyjaśniać te obserwacje eksperymentalne. Szczegółowo analizowany jest wpływ położenia, rozmiaru i momentu wprowadzenia źródła na zmiany w rozwoju somitów. Dodatkowo przebadano sytuację, w której w którymś obszarze zarodka pojawia się stały niedobór morfogenu, co otwiera nowe możliwości eksperymentalnego potwierdzenia poprawności sugerowanego mechanizmu somitogenezy. 

  k→l (N k + 1)P (N k + 1, N l -1) -d k→l (N k )P ] , (I.7) where d k→l is the probability of transition of a single molecule from cell number k to cell number l, where usually transition only between neighbouring cells is possible. By N k we denote population of molecules in cell number k. On the right-hand side of the above equation arguments of function P ({N }), that are the same as on the left-hand side, are omitted. Summation over j is done over all reactions, while over k and l over all spatial cells.

Fig. I. 1 (

 1 Fig. I.1(c).

  Fig. I.1(d)). According to linear stability analysis, there is no attraction or repulsion from the centre. Therefore nonlinear terms have to be taken into account to determine actual stability of the stationary state.

Figure I. 1 .

 1 Figure I.1.: Basic types of stationary states in a system of two variables x and y: (a) stable node, (b) saddle point, (c) stable focus, (d) centre. Eigenvectors of the Jacobian are denoted by v 1 and v 2 . (adapted from [26])

Figure

  Figure I.2.: An example function f (c) of a multistable system.

Figure I. 3 .

 3 Figure I.3.: An example probability distribution for the bistable Schlögl model. Dashed lines are stable state, dotted line is the unstable state.

  Fig. I.4.

Figure I. 4 .

 4 Figure I.4.: Supercritical Hopf bifurcation in the phase space of variables x, y. The bifurcation parameter is denoted by b. Bifurcation takes place for b = 0. Blue and red dots represent stable and unstable foci, respectively. (adapted from: http://www.me.rochester.edu/courses/ME406/webexamp6/bifurc6.pdf)

Figure I. 5 .

 5 Figure I.5.: Phase space of the Brusselator for parameter values A = 1, B = 3. Black line depicts the limit cycle, dotted lines are nullclines. The coloured points show results of stochastic simulations performed with the kinetic Monte Carlo method (cf. 4.2) for 1000 independent systems of volume V = 500 and initially in the state X = 3, Y = 0.93. Colours code position of the systems at different times: red at t = 11, blue at t = 22 and green at t = 33.

w = 8

 8 D/k. Example shapes of profiles of FKPP wave fronts are shown in Fig. I.6.

Figure I. 6 .

 6 Figure I.6.: Shape of a wave front in the FKPP model for two different propagation speeds v 1 > v 2 . The dashed line represents gradient of the profile at the inflection point. (adapted from [26])

  model (cf. Sec. 2.1.1) it is equal to

Figure I. 7 .Figure I. 8 .

 78 Figure I.7.: An example dispersion relation for the Gray-Scott model (dashed line).The vertical dotted lines depict wavenumbers that are allowed by boundary conditions.

r

  ∆t/2V collision pairs are randomly selected. The parameter v max r is the maximal relative velocity of a colliding pair and it can have any value as long as no pair of particles in the system has relative velocity greater than it. The formula above gives the number of collisions that would take place if relative velocities of all collision pairs would be equal to v max r . Collision of a given pair in fact takes place if v r > Rv max r

Figure

  Figure II.1.: A chick embryo with clearly visible somites (adapted from: http://www.maths.ox.ac.uk/groups/mathematicalbiology/research/pattern-formation).

Figure II. 2 .

 2 Figure II.2.: A schematic presentation of the "clock and wavefront" model with a 90 minutes cycle of oscillations. The blue colour denotes retinoic acid, greyfibroblast growth factor and in red are oscillations of gene expression in the PSM. Wave front propagation and the intermolecular clock are depicted on the left and right side of the notochord, respectively. In reality both this effects take place on the both sides of the notochord.

Równanie (I. 4 )

 4 opisuje układ jednorodny, tak więc uwzględnia jedynie fluktuacje całkowitych stężeń składników. Rozszerzając wyniki dla układów blisko równowagi często przyjmuje się, że amplituda fluktuacji jest odwrotnie proporcjonalna do pierwiastka z liczby cząsteczek, więc prawdopodobieństwo silnych fluktuacji, które spowodowałyby jakościowe zmiany w całym układzie, jest bardzo niskie w układach o rozmiarach makroskopowych. Takie globalne podejście do fluktuacji jest poprawne w przypadku układów liniowych, jednak zawodzi dla układów nieliniowych, w których lokalne fluktuacje mogą spowodować jakościowe zmiany w niewielkiej części objętości [5]. Następnie dyfuzja może spowodować rozprzestrzenienie się tej zmiany na całość układu. Prawdopodobieństwo takich dużych, lokalnych fluktuacji jest znacznie większe niż dużej fluktuacji w całej objętości ze względu na mniejszą ilość cząsteczek. Tak więc konieczne jest równanie master uwzględniające różnice w populacjach cząsteczek w różnych częściach systemu. Pozwala ono również badać układy, w których ze względu na warunki początkowe występują makroskopowe niejednorodności w stężeniach. Układ jest dzielony na komórki przestrzenne, między którymi zachodzi dyfuzyjna wymiana cząsteczek. Dla uproszczenia przyjmiemy, że w układzie znajduje się tylko jeden rodzaj cząsteczek. Interesuje nas prawdopodobieństwo P ({N }), gdzie {N } to zbiór populacji tych cząsteczek we wszystkich komórkach przestrzennych. Komórki muszą być na tyle małe, żeby można je było traktować jako jednorodne, pomimo makroskopowych niejednorodności w układzie. Z drugiej strony nie mogą być za małe, ponieważ transport pomiędzy komórkami jest opisywany przez dyfuzję, a nie przez ruch swobodny. Równanie master uwzględniające niejednorodności przestrzenne ma postać:

I. 7 )

 7 gdzie d k→l oznacza prawdopodobieństwo przejścia pojedynczej cząsteczki z komórki k do komórki l, przy czym uwzględnia się jedynie przejścia między sąsiednimi komórkami, a N k to ilość cząsteczek w k-tej komórce. Po prawej stronie powyższego równania pominięto te argumenty funkcji P ({N }), które są takie same jak po lewej stronie równania. Sumy po j brane są po wszystkich reakcjach, natomiast po k i po l są po wszystkich komórkach przestrzennych.

  Układy jednorodne są to układy, w których zaniedbywalne są przestrzenne różnice w stężeniach reagentów. Przykładem takiego układu są reaktory z mieszaniem. Zmienne zależą wówczas jedynie od czasu i nie uwzględnia się zjawisk transportu. Stan układu składającego się z N reagentów opisuje wektor c = [c 1 ,c 2 ,...,c N ], gdzie c i to stężenie i-tego reagenta. W ujęciu makroskopowym, dynamikę opisuje układ N równań różniczkowych:

Centrum -trajektorie wokół

  stanu stacjonarnego mają kształt elips (Rys. I.1(d)). Z analizy liniowej wynika, że nie ma przyciągania ani odpychania od centrum, więc w celu ustalenia faktycznej stabilności stanu potrzebne jest uwzględnienie członów nieliniowych. Rysunek I.1.: Podstawowe rodzaje stanów stacjonarnych w układach dwóch zmiennych x i y: (a) stabilny węzeł, (b) punkt siodłowy, (c) stabilne ognisko, (d) centrum . Wektory własne Jakobianu oznaczone są jako v 1 i v 2 . (za: Murray, Mathematical Biology [26]) 2.1.1. Multistabilność Multistabilność oznacza istnienie więcej niż jednego stabilnego stanu stacjonarnego. Sytuacja taka może mieć miejsce już w układach opisywanych jedną zmienną d t c = f (c). W tym przypadku warunkiem stabilności stanu stacjonarnego c s jest f ′ (c)| cs < 0. Przykładowa funkcja odpowiadająca układowi z trzema stabilnymi stanami stacjonarnymi pokazana jest na Rys. I.2. Stany niestabilne rozdzielają obszary atrakcji stanów stacjo-Rysunek I.2.: Przykładowa postać funkcji f (c) układu multistabilnego.

(I. 21 )

 21 Nadkrytyczna bifurkacja Hopfa następuje dla takiej zmiany parametrów, która powoduje przejście od ujemnej do dodatniej wartości śladu Jakobianu, przy zachowaniu det(J s ) > 0. Stan stacjonarny staje się niestabilny a układ opuści go w sposób oscylacyjny. Ponieważ nie ma już więcej stanów stacjonarnych, a żaden realistyczny schemat reakcji nie może pozwalać na dążenie stężeń do nieskończoności, układ osiągnie inny rodzaj atraktora. Okazuje się, że będzie to stabilny cykl graniczny otaczający badany stan stacjonarny. Charakterystyczne dla tego rodzaju bifurkacji jest to, że amplituda cyklu granicznego wzrasta wraz z oddalaniem się od bifurkacji. Schemat przejścia od stabilnego ogniska do stabilnego cyklu granicznego pokazany jest na Rys. I.4.Klasycznym modelem oscylacyjnego układu chemicznego jest "Brusselator"[32]. Za-Rysunek I.4.: Nadkrytyczna bifurkacja Hopfa w przestrzeni fazowej x, y. Jako b oznaczony jest parametr bifurkacyjny, bifurkacja zachodzi dla b = 0. Na niebiesko zaznaczone są stabilne ogniska, a na czerwono niestabilne ognisko. (za: http://www.me.rochester.edu/courses/ME406/webexamp6/bifurc6.pdf) wiera on dwa składniki o zmiennym stężeniu: X i Y.

  niczny jest stabilny wobec fluktuacji, poza przypadkiem fluktuacji równoległych do jego trajektorii, które nie będą zanikały. Efekt ten nazywa się dyfuzją fazy (phase diffusion) cyklu granicznego. Rozważmy duża liczbę identycznych układów, które początkowo są w tym samym punkcie na cyklu granicznym. Wraz z upływem czasu fazy ich oscylacji rozejdą się na skutek fluktuacji, aż dla t → ∞ ich fazy nie będą w ogóle skorelowane. To zjawisko zilustrowane jest na Rys. I.5. Już po jednym okresie oscylacji stany oscylatorów pokrywają cały cykl graniczny. Ponadto widać, że układ szybciej przemieszcza się przez obszar z dala od izoklin, natomiast wolniej blisko nich. Rysunek I.5.: Przestrzeń fazowa modelu Brusselatora dla parametrów A = 1, B = 3. Czarną linią ciągłą oznaczony jest cykl graniczny, a liniami kropkowanymi izokliny. Kolorowe punkty oznaczają wyniki symulacji stochastycznych wykonanych metodą kinetic Monte Carlo, przeprowadzonych dla 1000 niezależnych układów o objętości V = 500 i początkowo znajdujących się w stanie X = 3, Y = 0.93, który znajduje się na cyklu granicznym. Różne kolory oznaczają stany każdego z tych układów w różnych momentach: czerwone dla t = 11, niebieskie dla t = 22, a zielone dla t = 33.

33 )

 33 Szukamy rozwiązania tego równania, które by opisywało front przemieszczający się ze stałą prędkością v w kierunku dodatnim osi x [35]. Dodatkowo zakładamy, że front nie zmienia kształtu w czasie. Takie rozwiązanie ma postać: c(x,t) = c(x -vt) = c(z), z = x -vt. (I.34) Po podstawieniu do (I.33) otrzymujemy równanie różniczkowe zwyczajne drugiego rzędu: Dd zz c + vd z c + f (c) = 0. (I.35) Sposób rozwiązywania tego równania zależy od postaci funkcji f (c). Ważną grupę stanowią układy posiadające jeden stabilny stan stacjonarny dla c = 1 oraz jeden quasistabilny stan stacjonarny c = 0. Stan quasi-stabilny jest to taki stan stacjonarny, który z punktu widzenia dynamiki reakcji jest niestabilny, ale układ nie opuszcza go, ponieważ nie mogą pojawić się odpowiednie fluktuacje ze względu na brak cząsteczek odpowiedniego typu. Destabilizacja tego stanu może nastąpić tylko, gdy na skutek dyfuzji pojawią się takie cząsteczki. Front będzie się rozprzestrzeniał od obszaru w stanie c = 1 na obszar w stanie c = 0, co nakłada następujące warunki na funkcję c(z):

39 )

 39 Krytyczna szybkość frontu FKPP wynosi v = 2 √ kD a odpowiadająca jej szerokość w = 8 D/k. Przykładowe profile frontów FKPP pokazane są na Rys. I.6. Rysunek I.6.: Kształt front falowego w modelu FKPP dla dwóch różnych prędkości propagacji v 1 > v 2 . Linią przerywaną oznaczony jest gradient profilu frontu w punkcie przegięcia. (za: Murray, Mathematical Biology [26]) W przypadku, gdy pochodna f ′ (c) ma maksimum dla c > 0, szybkość frontu jest określana przez dynamikę reakcji nie na czele frontu, ale w jego wnętrzu. Takie fronty nazywamy frontami pchanymi. Nie ma ogólnego analitycznego wyrażenia na ich szybkość, ale również jest ona ograniczona od dołu przez wartość krytyczną v p , przy czym

  (patrz rozdział 2.1.1) prędkość ta wynosi

1 .

 1 0) > 0. (I.52) Biorąc pod uwagę (I.51) warunek ten może być spełniony tylko, gdy d = D v /D u = Tak więc struktury Turinga mogą powstawać tylko w układach, których składniki mają różne mobilności. Nie tracąc na ogólności, możemy przyjąć że d > 1. Kolejnym

  0, (I.62) co ogranicza możliwe struktury do tych o liczbie falowej k n = n π L , n = 1,2,.... Przyjmijmy parametry κ = 1, φ = 1, Y 0 = 5, dla których krytyczny stosunek współczynników dyfuzji wynosi d c ≈ 10.4. Załóżmy d = 15 i L = 40. Zależność dyspersyjna dla tych parametrów przedstawiona jest na Rys. I.7. Stan (X + ,Y + ) jest niestabilny wobec wszystkich modów z zakresu od k min do k max , przy czym warunki brzegowe nakładają dodatkowe ograniczenia na dostępne mody. Je-Rysunek I.8.: Rozwój struktury Turinga w jednowymiarowym układzie Graya-Scotta.Oś pozioma przedstawia wymiar przestrzenny, na osi poziomej znajduje się czas, natomiast kolor odpowiada stężeniu składnika X.żeli układ jest początkowo jednorodny, dla t → ∞ dominował będzie mod o największej wartości Reλ. Na Rys. I.8 przedstawiony jest rozwój struktury Turinga z układu jednorodnego.

6. 1 .

 1 Modelowanie powstawania struktur Turinga za pomocą techniki DSMCPodstawową trudnością przy mikroskopowym modelowaniu struktur Turinga jest uzyskanie wystarczająco różnych, co najmniej kilkukrotnie, stałych dyfuzji różnych składników. W mieszaninie sztywnych kul mobilność cząsteczek zależy od ich mas i rozmiarów.Aby dwa składniki miały różne dyfuzyjności, mieszanina musi składać się co najmniej z trzech rodzajów cząsteczek. Dwa spośród nich biorą udział w reakcjach, a trzeci pełni rolę rozpuszczalnika. Dyfuzyjności w takim układzie nie są stałe jak w analizie przedstawionej w rozdziale 2.2.2, lecz zależą od lokalnego składu mieszaniny[17]. Będą one jednak w przybliżeniu stałe w obszarze przestrzeni fazowej dookoła stanu, z którego powstaje struktura Turinga, jeżeli cząsteczek rozpuszczalnika jest znacząco więcej niż pozostałych składników. Wówczas dominują zderzenia pomiędzy reaktywnymi cząsteczkami a cząsteczkami rozpuszczalnika, natomiast zderzenia pomiędzy reaktywnymi cząsteczkami mają zaniedbywalny wpływ na ich dyfuzyjności. W symulacji, cząsteczki rozpuszczalnika pełnią także rolę rezerwuaru, z którego powstają cząsteczki reaktywne w stałym tempie, co podtrzymuje układ w stanie z dala od równowagi. Aby znacząco zwiększyć wydajność algorytmu, nie są symulowane zderzenia pomiędzy cząsteczkami rozpuszczalnika, które nie maja znaczenia z punktu widzenia struktury stężeń pozostałych składników. Symulacje przeprowadzone opracowaną metodą pokazały, że fluktuacje mogą powodować przejścia pomiędzy strukturami o różnych okresach.6.2. SomitogenezaZarodek zwierzęcy rozwija się z pojedynczej komórki, z której powstają kolejne, posiadające zdolność do różnicowania się w bardziej wyspecjalizowane typy komórek. Rozwój embrionu polega na skoordynowanym różnicowaniu się grup komórek, dzięki czemu powstaje skomplikowana struktura dojrzałego organizmu. W jaki sposób komórki uzyskują informację o swoim położeniu nie jest do końca wyjaśnione[76]. Stężenie niektórych substancji, zwanych morfogenami, ma wpływ na ścieżkę rozwoju niezróżnicowanych komórek. Dlatego też samoistne powstawanie uporządkowania przestrzennego w układach reakcja-dyfuzja, czyli struktury Turinga, uważa się za istotny mechanizm morfogenezy. Koronnym przykładem struktur Turinga w biologii są wzory na skórze zwierząt, jednak uważa się że mechanizm ten może być również odpowiedzialny za powstawanie innych struktur[26, 72, 77].Jedną z takich innych, stosunkowo prostych struktur, jest kręgosłup kręgowców. Na wczesnym etapie rozwoju struna grzbietowa zarodka dzieli się na regularne segmenty, zwane somitami (zob. Rys. II.1). Z somitów później rozwijają się kręgi i mięśnie kręgosłupa. W przybliżeniu, somitogeneza jest to proces polegający na powstawaniu jednowymiarowej struktury periodycznej. Sposób powstawania tej struktury nie został jeszcze wyjaśniony. Najczęściej postuluje się modele typu "clock and wavefront" (patrz rozdział 6.2.2)[78][79][80], jednak struktury Turinga też mogą wyjaśniać wiele charaktery-http://www.maths.ox.ac.uk/groups/mathematicalbiology/research/pattern-formation).

6. 2 . 2 .

 22 Schemat reakcja-dyfuzja dla modelu typu "clock and wavefront" Pierwszy model typu "clock and wavefront" został zaproponowany przez Cooke'a i Zeemana w celu wyjaśnienia somitogenezy[78]. Od tego czasu powstało wiele zbliżonych modeli, wszystkie jednak opierają się na założeniu, że struktura somitów powstaje na skutek interakcji pomiędzy oscylatorem wewnątrzkomórkowym (zegarem) a propagującym frontem chemicznym. Pierwotny model zakłada, że wewnątrz niezróżnicowanych komórek znajdujących się w tylnej części zarodka działają zsynchronizowane oscylatory chemiczne, natomiast od głowy w kierunku ogona zarodka przemieszcza się front chemiczny, który powoduje zatrzymanie się oscylatorów. W ten sposób powstaje struktura periodyczna, o okresie równym odległości przebytej przez front w trakcie jednego okresu oscylacji wewnętrznych zegarów.Rysunek II.2.: Schematyczne przedstawienie dwóch cykli mechanizmu "clock and wavefront". Niebieski kolor oznacza kwas retinowy, szary kolor oznacza czynnik wzrostu fibroblastów, natomiast na czerwono oznaczone są oscylacje w ekspresji genów w mezodermie przyosiowej. Po lewej stronie cewy nerwowej przedstawiona jest propagacja frontu, a po prawej -zegar wewnątrzkomórkowy. W rzeczywistości obydwa te efekty mają miejsce po obydwu stronach cewy nerwowej. Model ten zyskał popularność po odkryciu oscylacji ekspresji genów w tylnej części zarodków niektórych kręgowców. Matematycznie, model ten dotychczas był formułowany za pomocą makroskopowych równań nieposiadających podstaw mikroskopowych. W naszej publikacji przedstawiony jest pierwszy schemat reakcji elementarnych, który w połączeniu z dyfuzją składników daje zachowanie postulowane przez Cooke'a i Zeemana. Kryterium w trakcie opracowywania modelu była jego prostota -szukaliśmy minimalnego schematu reakcji, który by dawał pożądane zachowanie. Model składa się z dwóch powiązanych części -modelu antagonistycznego frontu falowego, opartego o model Schlögla oraz z modelu oscylatora chemicznego, zbliżonego do Brusselatora. Powiązanie pomiędzy nimi polega na tym, że składniki frontu chemicznego katalizują niektóre reakcje w modelu oscylatora. Dzięki temu po przejściu frontu chemicznego następuje jakościowa zmiana w dynamice oscylatora. Przed frontem, model oscylatora posiada atraktor w postaci stabilnego cyklu granicznego, natomiast po przejściu frontu ma on dwa stabilne stany stacjonarne. Dodatkowo, separatrysa między tymi stanami przecina cykl graniczny 50 II. Wyniki w warunkach przed przejściem frontu. Dzięki temu, po przejściu frontu komórki znajdą się w jednym z dwóch stanów stacjonarnych w zależności od tego, w jakiej fazie oscylatora znajdowały się w momencie nadejścia frontu chemicznego. W oparciu o opracowany model typu reakcja-dyfuzja przeprowadzone zostały symulacje metodą kinetic Monte Carlo, co pozwoliło dokładnie zbadać wpływ wewnętrznych fluktuacji na powstawanie struktury. Wykazały one, że dla powstania regularnej struktury niezbędna jest synchronizacja oscylatorów wewnątrzkomórkowych. W przypadku gdy synchronizacji nie ma, na skutek zjawiska dyfuzji fazy następuje desynchronizacja oscylatorów, a zatem otrzymana po przejściu frontu struktura jest nieregularna. Zaproponowano rozszerzony model, w którym sąsiednie oscylatory są ze sobą sprzężone za pomocą dyfuzji. III. Wnioski Przeprowadzona została szczegółowa analiza efektów stochastycznych występujących w niektórych ogólnych modelach układów nieliniowych. Zastosowano trzy poziomy opisu o różnej dokładności, przy czym każdy kolejny ujawniał nowe zachowania układu. W świetle przedstawionych wyników można wyciągnąć wniosek, że opis mezoskopowy, oparty na równaniu master, jest w większości przypadków wystarczająco dokładny. Opracowana została metoda pozwalająca zastosować opis mezoskopowy do układów, w których szybkie reakcje chemiczne powodują zaburzenie rozkładu prędkości cząstek. Zakres poprawności tej rozszerzonej metody został określony poprzez porównanie jej wyników z mikroskopowymi symulacjami komputerowymi, które wymagają długich obliczeń i mogą być trudne do zastosowania np. w biologii. Z drugiej strony makroskopowe podejście deterministyczne okazało się zbytnio uproszczone, ponieważ wewnętrzne fluktuacje mają znaczący wpływ na opisywane zjawiska. Opracowane zostały nowe techniki mikroskopowych symulacji komputerowych pozwalające modelować powstawanie struktur Turinga, które odgrywają istotną rolę w biologii rozwoju. Wyniki przeprowadzonych symulacji wskazują, że mechanizm Turinga, obok modelu "clock and wavefront", może tłumaczyć podział osi podłużnej zarodków kręgowców na segmenty. Zaproponowane zostały dodatkowe eksperymenty, które mogłyby rozstrzygnąć o poprawności któregoś z tych mechanizmów. Ponadto opracowany został model typu reakcja-dyfuzja dla mechanizmu "clock and wavefront", który stanowi pierwsze mikroskopowe ujęcie tego mechanizmu i może być pomocny przy budowie sztucznego układu o podobnym zachowaniu. Umożliwia on też badanie fluktuacji na poziomie mezoskopowym, co wskazało na konieczność sprzężenia oscylatorów w tym modelu. Warto również zaznaczyć, że najprostszy schemat reakcji konieczny do odwzorowania zachowania typu "clock and wavefront" jest znacząco bardziej skomplikowany od schematu, w którym mogą powstawać struktury Turinga.

  

  

  

  

  

  

  

  

  10) where J s is the Jacobian evaluated for the state c s . The solution of the set of Eqs. (I.10) is in the form of a linear combination of terms exp(λ i t), i = 1,...,N , where λ i are eigenvalues of matrix J s . Their respective eigenvectors are denoted by v i . A stationary state is stable if real parts of all eigenvalues are negative and unstable if at

least one is positive. If one of the eigenvalues has an imaginary part, the system in the vicinity of the state c s will move in the phase space in an oscillatory manner in the plane containing eigenvector of this eigenvalue. If all eigenvalues are real and have same sign, we call the stationary state a node. A stationary state that has only complex eigenvalues is called a focus. The most important types of stationary states in two-variable systems are shown in Fig. I.1. Depending on the eigenvalues of the Jacobian, they are

  Turing patterns, is considered as an important mechanism of morphogenesis. The most common example of Turing structures in biology are animal coat patterns, but it is believed that this mechanism can be also responsible for development of other patterns[26, 72, 77].

A rather simple example of another pattern is that of vertebrae. At an early stage of development, the anterior-posterior axis of an embryo is divided into regular segments, called somites (cf. Fig. II.1). Later, skeletal muscle, cartilage, tendons and vertebrae form from somites. In a simplified view, somitogenesis can be seen as formation of a onedimensional, periodic pattern. The mechanism of formation of this pattern is not clear.

Most often "clock and wavefront" type of models are proposed (cf. Sec. 6.2.2)

[78][79][80]

, but

  jakościowo różnych zachowań układu dynamicznego rośnie wraz z liczbą

	zywamy obszarem przyciągania tego atraktora. W dwuwymiarowej przestrzeni fazowej,
	jeżeli istnieje kilka atraktorów, ich obszary przyciągania rozdzielone są przez separatrysy
	(które przechodzą przez niestabilne stany stacjonarne) lub niestabilne cykle graniczne.
	Część przestrzeni fazowej, dla której jedna z pochodnych układu równań (I.1) równa zmiennych go opisujących. Moje badania dotyczyły efektów, które występują w układach jest zeru, nazywamy izokliną (nullcline). Stany stacjonarne znajdują się na przecięciach opisywanych jedną (multistabilność, fronty falowe) lub dwiema (chemiczne oscylacje, wszystkich izoklin. Stabilność stanu stacjonarnego można ustalić za pomocą liniowej struktury Turinga) zmiennymi, więc układy o większej liczbie zmiennych, w których analizy stabilności. Metoda ta oraz podstawowe rodzaje stanów stacjonarnych zostaną może wystąpić chaos deterministyczny [25], nie będą przedstawione. W kontekście chemicznym zmiennymi zazwyczaj są stężenia reagentów. W ujęciu ma-opisane w rozdziale. 2.1.
	kroskopowym dynamikę systemu opisuje układ równań deterministycznych na pochodne
	bilnym stanom stacjonarnym z równania makroskopowego, ale rzeczywisty stan układu
	nieustannie fluktuuje w okolicy tych maksimów. Stanom niestabilnym odpowiadają lo-
	kalne minima rozkładu prawdopodobieństwa.
	W ogólności punkt, lub zbiór punktów w przestrzeni fazowej, czyli przestrzeni zmien-
	nych opisujących układ, do którego dąży układ dynamiczny, nazywamy atraktorem.
	Układ o dowolnym stanie początkowym będzie ewoluował do któregoś z atraktorów dla
	t → ∞. Część przestrzeni fazowej, z której punkty zbiegają do danego atraktora, na-

poszczególnych stężeń po czasie, w postaci (I.1). Wartości stężeń, dla których wszystkie te pochodne równe są zeru, nazywamy stanami stacjonarnymi. Poza przypadkami krytycznymi, które odgrywają niewielką rolę w rzeczywistych układach, stan stacjonarny może być albo stabilny, tzn. gdy niewielkie wychylenia od niego zanikają, albo niestabilny, gdy niewielkie wychylenia powodują przejście do innego, stabilnego atraktora. W ujęciu deterministycznym układ nie jest w stanie opuścić stanu stacjonarnego.

W ujęciu mezoskopowym charakterystyczny dla układu jest stacjonarny rozkład prawdopodobieństwa, spełniający d t P ({N }) = 0. Niezależnie od początkowego rozkładu prawdopodobieństwa stanów, z biegiem czasu będzie on zbiegał do rozkład stacjonarnego. Rozkład stacjonarny posiada lokalne maksima dla wartości odpowiadających sta-

  ∂f 1 oznacza macierz Jacobiego liczoną w punkcie c s . Rozwiązanie układu równań (I.10) ma postać kombinacji liniowej czynników exp(λ i t), i = 1,...,N , gdzie λ i są to wartości własne macierzy J s , którym odpowiadają wektory własne v i . Stan stacjonarny jest stabilny, gdy część rzeczywista wszystkich tych wartości własnych jest ujemna, a nie-

						
	. . .	. . .	∂c N . . .	   	,	(I.10)
	∂f N ∂c 1	∂c N . . . ∂f N			
	gdzie J s					

stabilny, gdy któraś z nich jest dodatnia. Jeżeli któraś z wartości własnych ma niezerową część urojoną, układ w okolicy stanu c s będzie się poruszał po przestrzeni fazowej w sposób oscylacyjny w płaszczyźnie zawierającej wektor własny odpowiadający tej wartości własnej. Stan stacjonarny, którego wszystkie wartości własne są liczbami rzeczywistymi, nazywamy węzłem. Natomiast taki, którego wszystkie wartości własne są liczbami zespolonymi z niezerową częścią urojoną, nazywamy ogniskiem. Na Rys. I.1 pokazane są najważniejsze rodzaje stanów stacjonarnych w układach dwóch zmiennych. W zależności od wartości własnych Jakobianu, są to

  B, (I.13) gdzie kursywą oznaczone są stężenia danego składnika. Równanie to może mieć do trzech miejsc zerowych X 0 , X 1 , X 2 , odpowiadających dwóm stabilnym stanom stacjonarnym

przedzielonym przez jeden niestabilny stan stacjonarny. W ujęciu deterministycznym układ może przebywać w jednym ze stanów stacjonarnych i nie są możliwe przejścia między nimi. Bardziej interesujące zachowanie układu ujawnia się, gdy uwzględnimy efekty stochastyczne. Stacjonarny rozkład prawdopodobieństwa posiada lokalne maksima w punktach N X0 i N X2 , odpowiadające stabilnym stanom stacjonarnym. Pomiędzy tymi maksimami znajduje się lokalne minimum N X1 , czyli niestabilny stan stacjonarny.

  0→2 oznacza średni czas przejścia ze stanu N X0 do stanu N X2 , a τ 2→0 czas przejścia w odwrotnym kierunku. Jeżeli iloczyn znajdujący się na końcu równań (I.14) nie zawiera żadnych członów, uznaje się go za równy jedności.

	ξ) W układach opisywanych co najmniej dwiema zmiennymi możliwe jest pojawienie się , (I.15) nowego rodzaju atraktora nazywanego stabilnym cyklem granicznym. W przestrzeni fa-zowej jest to zamknięta orbita wzdłuż której podąża układ i do której zmierzają pobliskie trajektorie. W przebiegu czasowym oznacza to regularne i niegasnące oscylacje stężeń gdzie τ 2.1.2. Chemiczne oscylacje składników. Ogólne warunki na istnienie oraz liczbę cykli granicznych nie są znane. Na
	podstawie twierdzenia Poincarégo-Bendixsona wiadomo jednak, że w każdym jednospój-
	nym obszarze dwuwymiarowej przestrzeni fazowej, który nie zawiera stabilnych stanów
	stacjonarnych i do którego wnętrza wchodzą okoliczne trajektorie, musi istnieć cykl gra-
	niczny [30].

Przykładowy rozkład P (N X ) układu bistabilnego przedstawiony jest na Rys. I.3. Średni czas przejścia między stanami jest tym krótszy, im bliżej siebie znajdują się te stany oraz im większe jest prawdopodobieństwo odpowiadające stanowi niestabilnemu znajdującemu się między nimi. Rysunek I.3.: Rozkład prawdopodobieństwa stanów w układzie Schlögla. Linie kreskowane oznaczają stany stabilne, linia kropkowana stan niestabilny.

  Za pomocą analizy stabilności liniowej chcemy określić warunki, dla których stan ten jest stabilny wobec jednorodnych zaburzeń, ale niestabilny wobec zaburzeń przestrzennie periodycznych.

Analiza stabilności liniowej układu dwóch zmiennych bez dyfuzji opisana jest w rozdziale. 2.1.2. Na podstawie równania (I.

19

) wiemy, że stan (0,0) jest jednorodnie stabilny, jeśli wyznacznik macierzy Jacobiego w tym punkcie jest dodatni, a jej ślad jest ujemny:

  55) 

		0.2					
	Re(λ)	0.1					
		0					
		-0.1	0.6	k min	0.8	1	k max	1.2
					k		

  T . Prawdopodobieństwo, że przez czas τ nie zajdzie żaden proces, wynosi P (τ ) = exp (-ω T τ ). Jest to więc rozkład czasu oczekiwania na zajście jakiegokolwiek procesu. Długość kroku symulacji próbkowana jest z tego rozkładu i wynosi Gdy już zderzenie danej pary zostało zaakceptowane, konieczne jest określenie jego geometrii. Jako że mieszanina w komórce przestrzennej jest jednorodna, wszystkie parametry zderzenia są jednakowo prawdopodobne. Prędkości cząsteczek przed zderzeniem oznaczmy v 1 i v 2 , a ich prędkość względną jako v r = |v 1v 2 |. Przyjmijmy układ odnie-

	2, są wysoce wrażliwe na wewnętrzne fluktuacje. Prędkość ciągniętych frontów falowych zależy od dynamiki reakcji na czele W pierwszej części rozdziału II rozważany jest wpływ efektów stochastycznych na front FKPP i bistabilny układ Schlögla, w których szybkie reakcje chemiczne powodują zaburzenie rozkładu prędkości cząsteczek. Ważnym wynikiem tych prac jest zdefiniowa-nie równania master, które uwzględnia odchylenie rozkładu prędkości od Maxwellianu. Dokładność tego zaburzonego równania master określona jest poprzez porównanie jego 4. Wybrane metody numerycznego modelowania układów reakcja-dyfuzja Dotychczas zajmowaliśmy się analitycznym opisem układów chemicznych na różnych poziomach dokładności: deterministycznym makroskopowym i stochastycznym mezosko-czynając od przyjętych warunków początkowych, kolejne iteracje algorytmu dają stan układu w kolejnych, następujących po sobie momentach. Wymaga to przyjęcia kroku czasowego symulacji ∆t. Najprostszy jest algorytm Eulera, w którym układ ewoluuje według: c i (t + ∆t) = c i (t) + f i (c,t)∆t. (I.63) przybliżeniem gradientu pomiędzy dwiema komórkami jest ∇c x + ∆x 2 = c(x + ∆x) -c(x) ∆x . (I.69) Analogicznie, drugą pochodną przestrzenną w danej komórce liczymy jako Wybór metody modelowania zależy od gęstości układu, który chcemy badać. W ukła-dach gęstych zderzenia są częste oraz istnieją korelacje między położeniami i pędami par zderzeniowych. Dodatkowo, istotne mogą być dalekozasięgowe oddziaływania mię-jące się wewnątrz jednej komórki mogą zderzać się ze sobą. Taki podział pozwala również na uwzględnienie w symulacji niejednorodności występujących w układzie. Liczba zde-rzeń, pary molekuł które uległy zderzeniu oraz ich prędkości po zderzeniu losowane są we-4.3.2. Metoda dynamiki molekularnej Nazwa dynamika molekularna odnosi się do symulacji komputerowych ruchów wszyst-kich cząsteczek w mieszaninie, w których położenia i pędy zmieniane są zgodnie z pra-ich sumy ω ∆t = 1 ω T ln 1 r 2 , dzy cząsteczkami. W takim przypadku stosowana jest metoda dynamiki molekularnej dle rozkładów znanych z kinetycznej teorii gazu. Rozdział między ruchem a zderzeniami wami dynamiki Newtona [55]. Jest to więc metoda w pełni deterministyczna. Ważną (I.71) (MD), która polega na stosowaniu równań Newtona do położeń i pędów wszystkich czą-może być stosowany jedynie, gdy krok czasowy ∆t jest znacznie krótszy od średniego sienia środka masy pary zderzeniowej, poruszający się z prędkością v m . W tym układzie rolę odgrywa dobór potencjałów międzycząsteczkowych U (r). Ogół sił niezwiązanych z frontu, gdzie cząsteczek jest mało, a więc fluktuacje mają duże znaczenie dla rozwoju całego frontu. W układach bistabilnych spontaniczne przejścia między stanami stacjonar-nymi wywoływane są przez fluktuacje. W przypadku oscylatorów chemicznych, fluktu-acje powodują dyfuzję fazy. Z kolei struktury Turinga mogą rozwinąć się z przestrzennie jednorodnego stanu stacjonarnego, gdy stochastyczne zaburzenie wyprowadzi układ z tego stanu. Kolejną przyczyną, dla której ważne jest uwzględnienie fluktuacji w opisie, jest to, że wymienione efekty odgrywają istotną rolę w funkcjonowaniu komórek orga-nizmów żywych [26]. Fronty chemiczne pozwalają na rozprzestrzenianie się informacji wewnątrz i pomiędzy komórkami, cykle graniczne stanowią zegary komórkowe, komórki mogą przebywać w jednym z wielu stanów stacjonarnych, natomiast struktury Turinga są jednym ze sposobów powstawania zróżnicowania przestrzennego komórek. Ponieważ populacje cząsteczek wewnątrz komórek są niewielkie, w analizie powyższych efektów Celem pracy jest zbadanie wpływu fluktuacji wewnętrznych na zachowanie wybranych układów nieliniowych przy pomocy analizy na poziomach mezoskopowym i mikrosko-powym. Punktem odniesienia będzie makroskopowe podejście deterministyczne. Poziom skomplikowania układów dynamicznych jest na tyle duży, że do ich analizowania trzeba użyć metod numerycznych, których opracowanie stanowi istotną część tej pracy. Wydaj-nym sposobem symulacji uwzględniających wewnętrzne fluktuacje jest technika kinetic Monte Carlo, oparta na chemicznym równaniu master. Jednakże metoda ta zawiera sze-reg przybliżeń i konieczne jest potwierdzenie jej wyników przez symulacje na poziomie mikroskopowym. Oddzielną grupą technik komputerowych są symulacje mikroskopowe. Nie wywodzą się one z metod analitycznych, lecz mają na celu możliwie najdokładniejsze odwzorowa-nie ruchów pojedynczych cząsteczek. Są to metody najbliższe rzeczywistym układom, ale też najbardziej wymagające obliczeniowo. Dwie metody mikroskopowe, dynamika molekularna i direct simulation Monte Carlo, przedstawione są w rozdziale 4.3. 4.1. Numeryczne rozwiązywanie równań makroskopowych Istnieje wiele metod iteracyjnego rozwiązywania układów równań różniczkowych zwy-policzenia prawdopodobieństw przejścia dla wszystkich M procesów ω i , i = 1,...,M oraz cząsteczek. Dlatego wprowadza się podział układu na komórki i tylko cząsteczki znajdu-układów i krótkich czasów symulacji. czajnych, takich jak równania makroskopowe dla układu jednorodnego (I.9) [49]. Za-O (∆t 4 ). W przypadku układów reakcja-dyfuzja, które opisują równania różniczkowe cząstkowe w postaci (I.32), również można użyć powyższych metod, jednak najpierw konieczne jest przybliżenie tych równań za pomocą równań różniczkowych zwyczajnych [51]. Jednym ze sposobów na to jest metoda różnic skończonych. Opiera się ona na dyskretyzacji przestrzeni na komórki o rozmiarach ∆x. Dla układu jednowymiarowego najprostszym komórki musi być starannie dobrany. Powinna być ona na tyle mała, żeby można ją było traktować jako jednorodną, ale na tyle duża, żeby transport między komórkami miał charakter dyfuzyjny. Ponadto wraz ze zmniejszaniem rozmiaru komórki znacząco spada efektywność symulacji. Zainicjowanie symulacji polega na zdefiniowaniu liczby cząsteczek w poszczególnych komórkach. Krok symulacji polega na zajściu pojedynczego procesu w układzie. W ukła-dach reakcja-dyfuzja takim procesem może być pojedyncza reakcja elementarna, bądź dyfuzyjny ruch jednej cząsteczki do sąsiedniej komórki. Krok symulacji zaczyna się od czających i zderzających się ze sobą cząsteczek, dlatego też najdokładniejsze metody modelowania komputerowego opierają się na śledzeniu położeń i pędów pojedynczych cząsteczek, co nazywamy modelowaniem mikroskopowym [54]. Takie podejście jest naj-bliższe rzeczywistości i pozwala na lepsze zrozumienie zjawisk trudnych do uchwycenia eksperymentalnie. Stosowanie bardziej uproszczonych metod mezoskopowych czy ma-kroskopowych może być uzasadnione przez porównanie z modelowaniem mikroskopo-wym. Podstawową wadą algorytmów mikroskopowych jest wymagany bardzo duży na-kład pracy komputerowej, dlatego ich stosowanie jest ograniczone tylko do niewielkich Direct simulation Monte Carlo jest to metoda symulacji mikroskopowych dla układów o niewielkich gęstościach opracowana przez Birda [18,57]. Może być ona traktowana jako rozwiązywanie nieliniowego równania Boltzmanna [58]. Została stworzona do modelowa-nia przepływów w trakcie lotów w wysokich partiach atmosfery, jednak znalazła również zastosowanie w takich dziedzinach jak nanoprzepływy i badanie fluktuacji chemicznych. Opiera się na rozdzieleniu ruchu cząsteczek od zderzeń między nimi. W każdym kroku symulacji ∆t położenia zmieniane są zgodnie z równaniami ruchu r(t+∆t) = r(t)+v∆t. Dana cząsteczka, przebywając odległość v∆t, mogła zderzyć się z którąś z okolicznych Parametr v max r wszystkich par wynosiły v max r runkiem v r > Rv max r cząsteczek w komórce. każdej komórki. narnych i jest bardziej efektywny od sprawdzania zderzeń pomiędzy wszystkimi parami od zderzeń i ich liczba w kroku czasowym jest losowana z rozkładu dwumianowego dla [0,1). Taki sposób losowania zderzeń jest zgodny z Boltzmannowską całką zderzeń bi-czającej energii faktycznie będzie reaktywne. Reakcje jednocząsteczkowe są niezależne , gdzie R to liczba losowa o rozkładzie jednorodnym na zakresie ryczny. Jest to prawdopodobieństwo z zakresu (0,1] na to, że dane zderzenie o wystar-. Zderzenie wybranej pary rzeczywiście zachodzi pod wa-wewnętrzną strukturą cząsteczek a zajściem reakcji, wprowadza się tzw. czynnik ste-ści względnej. Powyższe wyrażenie określa liczbę zderzeń, gdyby prędkości względne dla jest większa od energii aktywacji reakcji. Aby uwzględnić ewentualną zależność między być wybrany dowolnie, o ile żadna para cząsteczek nie ma od niego większej prędko-czym takie zderzenie uznaje się za elastyczne. Reakcja zachodzi, gdy energia zderzenia oznacza maksymalna względną prędkość między cząsteczkami i może Zderzeniowe reakcje chemiczne polegają na zmianie typu chemicznego cząsteczek, przy Algorytm Rungego-Kutty jest numerycznie stabilny, a jego całkowity błąd jest rzędu pulacji poszczególnych składników we wszystkich komórkach przestrzennych. Rozmiar Na fundamentalnym poziomie układ reakcja-dyfuzja jest to mieszanina przemiesz-4.3.1. Metoda Direct Simulation Monte Carlo nieliniowych w kontekście biologicznym powinny być uwzględniane fluktuacje. wyników z symulacjami mikroskopowymi. Druga część rozdziału II opisuje badania nad wpływem fluktuacji na powstawanie struktur przestrzennych, z naciskiem na powstawanie kształtu w trakcie rozwoju embrio-nalnego. Opracowana została metoda symulacji mikroskopowych powstawania struktur Turinga. Następnie była ona stosowana w celu objaśnienia wyników eksperymentów, polegających na zaburzaniu rozwoju zarodka. W ten sposób pokazano, że mechanizm Turinga jest zgodny z tymi doświadczeniami. Ponadto opracowano prosty schemat re-akcji dla modelu "clock and wavefront", który wyjaśnia powstawanie jednowymiarowych struktur periodycznych w zarodkach kręgowców. Na podstawie tego schematu możliwe było zastosowanie metody kinetic Monte Carlo, aby zbadać wpływ fluktuacji na model "clock and wavefront". powym. W obydwu przypadkach, poza najprostszymi modelami, niemożliwe jest anali-tyczne rozwiązanie równań dynamiki układu. Aby uzyskać informacje o ewolucji układu, niezbędne jest więc zastosowanie metod numerycznych. Najbardziej efektywne obliczeniowo jest rozwiązywanie układów równań makroskopo-wych. Algorytmy do tego służące opisane są w rozdziale 4.1. Za ich pomocą uzyskujemy deterministyczną trajektorię układu, zależną jedynie od warunków początkowych. W ujęciu mezoskopowym, bezpośrednie całkowanie numeryczne równania master jest zazwyczaj nieefektywne ze względu na ogromną liczbę stopni swobody. Dlatego do jego rozwiązywania stosuje się algorytm kinetic Monte Carlo (KMC), opisany w rozdziale 4.2. Takie symulacje dają pojedynczą stochastyczną trajektorię systemu, co odpowiada prób-kowaniu z rozkładu prawdopodobieństwa stanów. Dlatego uzyskanie pełnej informacji o tym rozkładzie wymaga powtórzenia wielu tego rodzaju symulacji. Metoda ta posiada szereg wad: jest numerycznie niestabilna i jej błąd całkowity jest proporcjonalny do ∆t. Bardziej dokładna, a co za tym idzie pozwalająca na obranie większych kroków czasowych, jest metoda Rungego-Kutty [50]. Dla układu jednoskład-nikowego polega ona na liczeniu stanu układu w kolejnych krokach jako: c(t + ∆t) = c(t) + 1 6 ∆t(k 1 + 2k 2 + 2k 3 + k 4 ) (I.64) k 1 = f (c,t) (I.65) k 2 = f c + ∆t 2 k 1 ,t + ∆t 2 (I.66) k 3 = f c + ∆t 2 k 2 ,t + ∆t 2 (I.67) k 4 = f (c + ∆tk 3 ,t + ∆t) . (I.68) ∆c (x) = c(x -∆x) + c(x + ∆x) -2c(x) ∆x 2 . (I.70) 4.2. Modelowanie mezoskopowe -metoda kinetic Monte Carlo W oparciu o równanie master (patrz rozdział 1.2) Gillespie stworzył algorytm służący do modelowania stochastycznej trajektorii układu, nazywany kinetic Monte Carlo [52, 53]. Zaletą tej mezoskopowej techniki jest to, że przy stosunkowo niskim nakładzie pracy komputera, uwzględnia wewnętrzne fluktuacje układu bez żadnych dodatkowych przybli-żeń dotyczących ich rozkładu. Z drugiej strony ma ona wszystkie ograniczenia równania jest znajomość wyrażenia na prawdopodobieństwa przejść występujące w tym równaniu. W przypadku układu niejednorodnego przestrzennie stan układu opisuje wektor po-4.3. Modelowanie mikroskopowe pojedynczych reakcji nie są istotne, dlatego w symulacjach używamy znacznie uprosz-v * 2 = v m -r 2 . (I.75) v * czonych modeli reakcji, a same cząsteczki są traktowane jako kule. v * 1 = v m + r 2 (I.74) v * master, czyli może być użyta tylko do symulowania procesów Markowa oraz potrzebna steczek [55]. Metoda ta jest nieefektywna dla układów rozrzedzonych, w których średnia czasu swobodnego. Wtedy prawdopodobieństwo kilku zderzeń pojedynczej cząsteczki w elastyczne rozpraszanie w zderzeniu twardych kul jest izotropowe i rozkład prawdopo-wiązaniami chemicznymi ani oddziaływaniem elektrostatycznym między jonami nazy-gdzie r 2 to liczba losowa z rozkładem jednorodnym na przedziale (0,1]. Prawdopodo-droga swobodna jest znacznie większa od rozmiarów cząsteczek. W takich przypadkach czasie jednego kroku jest znikomo małe. dobieństwa kierunku prędkości względnej po zderzeniu v * r jest sferycznie jednorodny. wamy oddziaływaniem Van der Waalsa. Składa się na nie krótkozasięgowe odpychanie bieństwo, że po czasie ∆t zajdzie i-ty proces wynosi ω i /ω T . W celu wybrania procesu bardziej adekwatna jest metoda direct simulation Monte Carlo (DSMC) [18]. Średnią liczbę zderzeń wewnątrz komórki przestrzennej można policzyć jako iloczyn Współrzędne wektora względnej prędkości cząsteczek po zderzeniu wynoszą [18] wynikające z reguły Pauliego oraz dalekozasięgowe przyciąganie między stałymi lub in-losowana jest liczba r 1 z rozkładu jednorodnego na przedziale (0,1]. W danym kroku czasowym zajdzie proces numer β, dla którego spełnione jest β-1 i=1 ω i < r 1 ω T ≤ β i=1 ω i . (I.72) Na koniec kroku symulacji zmieniany jest wektor stanu układu zgodnie z wylosowanym procesem oraz zwiększany jest czas symulacji o ∆t. Istotną kwestią jest sposób modelowania reakcji chemicznych zachodzących w ukła-dzie. Ideałem byłoby uwzględnianie orbitali każdej z cząsteczek i przeprowadzanie re-akcji zgodnie z zasadami mechaniki kwantowej [56]. Jednak takie szczegółowe podejście jest wymagające obliczeniowo nawet w przypadku tylko kilku cząsteczek, podczas gdy przedstawiona praca poświęcona jest badaniu kolektywnego zachowania układów zawie-rających dziesiątki tysięcy molekuł. Z takiego globalnego punktu widzenia szczegóły dukowanymi dipolami. Najczęściej stosowanym przybliżeniem tych oddziaływań między liczby par cząsteczek w komórce N (N -1)/2, gdzie N to liczba cząsteczek, i prawdopo-dobieństwa zderzenia danej pary. To prawdopodobieństwo zależy od przyjętego modelu oddziaływań. Tutaj ograniczymy się do modelu twardych kul, jednak w zastosowaniach inżynieryjnych używa się bardziej dokładnych modeli [18, 59]. Zderzenie następuje, gdy v * dwiema cząsteczkami jest potencjał Lennarda-Jonesa w postaci [60]: r = (cos(β)v sienia systemu liczymy jako U LJ (r) = 4ǫ σ r 12 -σ r 6 , (I.76)

odległość między środkami cząsteczek z pary zderzeniowej będzie mniejsza od sumy ich promieni, oznaczonej jako σ. W układzie odniesienia jednej z cząsteczek z pary, druga z nich w kroku czasowym przeczesuje objętość V coll = πσ 2 v r ∆t, gdzie v r to prędkość względna między nimi. Ponieważ zakładamy, że komórka przestrzenna jest jednorodna, prawdopodobieństwo zderzenia między parą cząsteczek wynosi V coll /V , gdzie V to objętość komórki. W każdej komórce losuje się N (N -1)πσ 2 v max r ∆t/2V par zderzeniowych. r , sin(β) cos(ǫ)v r , sin(β) sin(ǫ)v r ) , (I.73) gdzie cos(β) jest losowane z rozkładu jednorodnego na przedziale [-1,1], a ǫ z rozkładu jednorodnego na przedziale [0,2π). Prędkości cząsteczek po zderzeniu w układzie odniegdzie r to odległość między środkami cząsteczek, ǫ to głębokość minimum potencjału, a σ to odległość, w której potencjał równy jest zeru. Przyciąganie międzycząsteczkowe jest stosunkowo słabe i w niektórych przypadkach może być pominięte, np. w wysokiej temperaturze lub w niezbyt gęstym układzie. W takim przypadku, gdy interesuje nas wyłącznie odpychanie, można używać potencjału Weeksa-Chandlera-Andersena [61]:

  1/6 σ = r ir j . Wypadkowa siła działająca na j-tą cząsteczkę wynosi F j = i =j F ji .W każdym kroku symulacji liczona musi być wypadkowa siła działająca na każdą cząsteczkę w układzie. Jest to część algorytmu wymagająca największego nakładu pracy komputera, który jest proporcjonalny do kwadratu całkowitej liczby cząsteczek. Aby uzyskać skalowanie liniowe wraz z wielkością układu stosuje się przybliżenie, polegające na ucięciu potencjału międzycząsteczkowego w odległości, dla której jego gradient jest zaniedbywalnie mały. Dla potencjału Lennard-Jonesa jest to zazwyczaj r c = 2.5σ. Symulacje z zastosowaniem potencjału WCA są szybsze, ponieważ dla niego r c = 2 1/6 σ, więc sąsiedztwo każdej cząsteczki jest mniejsze. Po policzeniu sił należy numerycznie scałkować równania ruchu. Niektóre metody do tego służące opisane są w rozdziale 4.1. W

	gdzie r			
	Jest to potencjał LJ ucięty w odległości odpowiadającej jego minimum. Dodatkowo
	jest on przesunięty o ǫ, aby zachować ciągłość potencjału.	
	Siła działająca na cząsteczkę i na skutek oddziaływania z cząsteczką j wynosi	
	F ij = -d r U (|r|)	r |r|	,	(I.77)

najprostszym przypadku zmieniamy pędy cząsteczek o ∆v = F ∆t/m, a ich położenia o ∆x = v∆t. Długość jednego kroku ∆t jest zazwyczaj rzędu femtosekund. Jeżeli chcemy symulować układ kanoniczny, czyli pozostający w kontakcie termicznym z otoczeniem, konieczne jest wprowadzenie algorytmu symulującego termostat. Jednym z popularnych sposobów jest termostat Nosé-Hoovera [62]. Przestrzeń fazową układu rozszerza się o dodatkową zmienną η. Jej równanie ruchu ma postać d t η = E k -gk B T, (I.78) gdzie E k to całkowita energia kinetyczna cząsteczek, g to liczba stopni swobody, a T to temperatura otoczenia. Do siły działającej na każdą cząsteczkę dodaje się człon odpowiadający wymianie energii z otoczeniem, tak że całkowita siła działająca na j-tą cząsteczkę równa jest

  Zaburzenie powoduje przesunięcie obszaru bistabilności w przestrzeni parametrów modelu. Oznacza to, że układ, który przy tradycyjnej analizie nieuwzględniającej zaburzenia byłby monostabilny, w rzeczywistości jest bistabilny. Możliwe są też przejścia od bistabilności do monostabilności, oraz pomiędzy układami monostabilnymi, ale o znacząco różnych położeniach stanów stabilnych. Wykazano też możliwość zbliżenia się stanów stabilnych do siebie, co powoduje zmniejszenie średniego czasu przejścia między nimi o kilka rzędów wielkości. Wyniki uzyskane z równania master uwzględniającego zaburzenie, które zostało wyprowadzone w tej pracy na podstawie równania makroskopowego[68], zostały porównane z symulacjami metodą DSMC. Zaburzone równanie master przestaje być poprawne dla energii aktywacji mniejszych niż 6k B T .

	5. Wpływ zaburzenia rozkładu prędkości cząsteczek dział 2.1.1).
	na dynamikę układu
	Szybkość reakcji chemicznych można obliczyć na podstawie częstotliwości zderzeń mię-
	dzycząsteczkowych, przy czym zazwyczaj przyjmuje się Maxwellowski rozkład prędkości
	cząsteczek poszczególnych reagentów [17, 63]. W większości przypadków takie założenie
	jest zasadne, jednak gdy reakcja jest bardzo szybka, t.j. gdy częstotliwość zderzeń reak-
	tywnych jest porównywalna z częstotliwością zderzeń elastycznych, pojawia się znaczące
	odchylenie rozkładu prędkości reagentów od rozkładu równowagowego [64-68]. Spowo-
	dowane jest to z tym, że cząsteczki o wyższej energii mają większe prawdopodobieństwo
	przereagowania od tych o niższej energii. Efekt ten ma dwie przyczyny: wzrost często-
	tliwości zderzeń wraz ze wzrostem energii kinetycznej oraz istnienie energii aktywacji
	Zaburzenie rozkładów prędkości ma wpływ na dynamikę reakcji, która odbiega od
	standardowego równania na szybkość reakcji. Takie równanie zawierające poprawki może
	być wyprowadzone z równania Boltzmanna za pomocą rachunku zaburzeń przy zachowa-
	niu większej liczby członów rozwinięcia [69]. W ogólności uzyskujemy poprawki zarówno
	dla szybkości reakcji jak i dyfuzji składników.
	5.1. Multistabilność
	Zaburzenie rozkładu prędkości może wpływać na liczbę stanów stacjonarnych układu,
	co pokazano na modelu Schlögla, który jest typowym modelem bistabilności (patrz roz-

reakcji, t.j. minimalnej energii pary zderzających się cząsteczek, poniżej której nie mogą one przereagować. W efekcie obserwuje się przesunięcie rozkładu prędkości produktów w stronę wysokich energii a substratów w stronę niskich energii. Jeżeli układ jest termoneutralny, t.j. nie ma efektu cieplnego reakcji, sumaryczny rozkład prędkości wszystkich cząsteczek w układzie nie jest zaburzony i ma postać rozkładu równowagowego.

Acknowledgements iv I. Introduction Acknowledgements Most of all, I would like to express gratitude for my supervisors, Bogdan Nowakowski and Annie Lemarchand, for sharing their extensive knowledge, encouraging my research and for their understanding. During my studies I spent an inspiring period of time at the Université Pierre et Marie Curie, mostly thanks to Annie Lemarchand, Charles Antoine and Laurence Signon, thank you. I was equally well greeted at the Roskilde University, for which I am thankful to members of the Glass and Time group, especially Jesper Schmidt Hansen and Claire Lemarchand.

I am indebted to my colleagues from the Institute of Physical Chemistry, Agnieszka

This dissertation was prepared within the International PhD in Chemistry Studies at the Institute of Physical Chemistry of the Polish Academy of Sciences in Warsaw. Programme co-financed by the European Regional Development Fund, Operational Program Innovative Economy 2007-2013. The dissertation was prepared within the co-tutorial agreement between the Institute of Physical Chemistry of the Polish Academy of Sciences in Warsaw and the Université Finally, I would like to thank Ola for her constant support and never doubting my choices. The Project was operated within the Foundation for Polish Science International PhD Projects Programme co-financed by the European Regional Development Fund, Operational Program Innovative Economy 2007-2013. Additional financial support during my stays in France was provided by the "doctorat en cotutelle" programme of the French Government. iv

Spis treści

Streszczenie i

Podziękowania v

I. Wprowadzenie

Każdej liczbie falowej spełniającej h(k 2 i ) < 0 odpowiada przynajmniej jedna dodatnia wartość własna λ. Zależność większej z dwóch Reλ od liczby falowej k i nazywa się zależnością dyspersyjną (dispersion relation).

Jednym z najprostszych modeli, w których mogą powstać struktury Turinga, jest model Graya-Scotta [46,47]. Uwzględnia on dwie reakcje zachodzące w reaktorze przepływowym

Łącznie reakcje te odpowiadają zamianie substratu Y w produkt P . W reaktorze wywołany jest przepływ φ, który dostarcza składnik Y w tempie φY 0 i usuwa składniki X, Y i P. Równania reakcja-dyfuzja mają postać

Skala czasu została dobrana tak, aby uzyskać κ 1 = 1, wtedy κ 2 = κ, natomiast skala odległości tak, aby D X = 1, wtedy D Y = d. Układ może mieć trzy stany stacjonarne (0,Y 0 ),(X -,Y -),(X + ,Y + ), gdzie

Rozważmy powstawanie struktury Turinga ze stanu (X + ,Y + ) w jednowymiarowym układzie Graya-Scotta o długości L. Jakobian w tym stanie wynosi