Expansion of the master equation for a biomolecular selection model

Bullettn ofMathemotrcalBio/ogy

Vol. 47. No. 2, pp. 23 l-238.

1985.

Printed in Great Britain

OU92-8240/X5$3.00

+ 0.00

Pergamon Press Ltd. Q 1985 Society for Mathematical

Biology

EXPANSION OF THE MASTER EQUATION FOR A BIOMOLECULAR SELECTION MODEL n

H. K. LEUNG Physics Department, National Central University, Chung-ii, Taiwan 320, Republic of China

A stochastic model based on Eigen and Schuster’s theory of biomolecular self-replication is studied by treating the master equation with the system-size expansion technique. The steady-state results are found to be in good agreement with the previous results and with those. derived from the principle of detailed balancing. Multispecies competition and coexistence are studied carefully with the conclusions that a stable steady state is predicted for the former and a metastable state for the latter. The stochastic selection processes are also analyzed and discussed.

Recently, Eigen and Schuster’s theory of biomolecular 1. Introduction. self-replication (Eigen, 197 1; Eigen and Schuster, 1979) was studied stochastically (Jones and Leung, 1981; Leung, 1984). Owing to the nonlinear nature of the model, the corresponding master equation is practically unsolvable. Even the moment equations have to be truncated so that the average molecular numbers and the variances can be studied approximately. Most stochastic behaviors are indeed included in the first two moments, however, it would be more desirable to evaluate the probability distribution function from which more stochastic aspects of the self-replication processes will be revealed. An alternate approximation is needed to achieve this goal. We note also that the moment-expansion scheme (Jones and Leung, 198 1) involves the neglect of the third correlation. This and higherorder correlations have been recently found to be important, for some cases, even in the linear system (Schuster and Sigmund, private communication). In this paper, results from different approximation schemes will be compared. The model system under consideration consists of N species of biomolecular information carriers which are capable of replicating themselves accurately. External constraint is imposed so that the carrying capacity of the system is Jz, while the total number of molecules in the system is variable. The deterministic equations describing the self-replication processes are given by,

dnildt = g(ni) - I),

(1)

with

id%) = At+%, rj({n}) = Djnj

+r:wfnp/n.

(3) 231

232

H.K.LEUNG

In the above, ni is the number of molecules for the ith species and Wi = Aj - Q is the selective value. Aj and Di are respectively the self-replication

and the degradation rates. In equation (I), the complicated processes of molecular replication are treated as the single-step autocatalytic reaction, so that the master equation approach developed for chemical dynamics (McQuarrie, 1967) can then be employed to describe the stochastic behaviors of the model system. The probability of having nj molecules of species j at the time t is described by the probability function P({n}; t) which satisfies the master equation (Gardiner, 1983), aP({n}; t)/& =

&(nj -

1>P(n, . . . nj -

1 . . . ,,N; t)

i

+ rj(n

1.

.

.

nj + 1 . . . nN)P(nl . . . n/ + 1 . . . nN; t) -

[g(q) + g{n)>lP({n); t>).

(3)

This stochastic description is valid for a system with large S2 since a correct stochastic theory must reduce to the deterministic formulation in the thermodynamic limit of C! + 00 (Oppenheim et al., 1969; Kurtz, 1972). We have shown previously, with the moment expansion approximation scheme, that this master equation (3) produces approximately the deterministic rate equation together with fluctuations which are insignificant if a is large (Jones and Leung, 198 1; Leung, 1984). These facts also enable us to employ the system-size expansion approximation (Van Kampen, 1976; Van Kampen, 1981) to expand the master equation (3). This new treatment allows the determination of the otherwise mathematically inaccessible but biologically interesting regime, and provides new results which can be used to compare with those derived from a different approximation scheme. The system-size expansion tech2. Expansion of the Master Equation. nique was designed to treat the stochastic problems with small relative fluctuations (Van Kampen, 1976). Each variable n,(t) can be split into deterministic and stochastic parts, q(t)

= L?r#q(t) + W2x,(t).

(4)

The probability function P((n}; t) now goes over to II({x);t) which describes the probability of having the fluctuating component xl(t) while its deterministic part is &l(t). These two probability functions are related by, I-I({x}; t) = swP({n}; After substitution

t).

and expansion, equation (3) reduces to

(5)

ABIOMOLECULARSELECTIONMODEL

233

(6)

where cu(k)is the reduced jump moment of order k for the species j, c4k) = [g(nj) + (-l)krj({n})l 1

/a.

(7)

Wj#j4i

(8)

From equation (6), we find for the leading terms, d@i/d t = c~i(l)= Wi$i -x

I

which is just the deterministic rate equation = C&i(t). The next order terms give,

for the average population

iii(t)

(9) which is the multivariate Fokker-Planck equation. The solution of equation (9) can be found through lengthy mathematics to be a multivariate Gaussian function with parameters too long to be listed here. We shall illustrate it in the next section for a simple case. Here it would be helpful to put down the stochastic equations for the moments of xi(t). The averages and the variances Uii = Z$j - ZiX/ are found to satisfy the following, dxi/d

dikxk ,

t =x

(10)

k

d oil/d t = x (dik okj +d_jk oki) +ai&ij

(11)

k

where,

dij The real fluctuations by OijS2.

= ad) I /a#1,

and correlations

@i =

01”’

(12)

in the molecular numbers nj are given

As we shall see later the simple 3. Constrained One-species Self-replication. one-species system shares the same stochastic features of the steady state (SS) common to the multispecies systems. For N = 1, the Fokker-Planck equation (9) can be simplified as,

234

H.K.LEUNG

an/at

= --w(i - 2#)a(dyax

where p = A/W is the metabolism solution to the equation (8),

+ w/2(2jd - i + 2#)a2n/ax2,

(13)

factor, and 4 = ii/s2 is the deterministic

N) = &A& + ( 1

-

40)

ew

(- W>-‘,

(14)

with #o = #J(O). The stochastic importance of p is evident from the FokkerPlanck equation (13) that 1 enters into the diffusion term which plays a role in diffusing the.II(x; t) curve while the deterministic type drift term tries to have II(x; t) peak at x = 3. The probability function can then be solved from equations (13) and (5) and is given by the normalized Gaussian form, P(n; t) =~(2do(t))

exp {--In - Ch#4t)l2/2s2a(t)}.

(15)

In the above, the mean variance u(t) is a solution of equations (10) and (11): o~t~=92~1-~~2{~~6~-11)~/~~-~~-~6~--~~~~l~~-~~/91 -

c&4 -

1x1 -

@>I# + P92/(1 -

~)‘I$,>.

(16)

The ‘process of the one-species replication is now described by the Gaussian P(n; t) which peaks at the average population ii(t) = C24(t) and spreads with a narrow width equalling the fluctuation 4 [ Sh(t)] . At any time there is a finite, though very small, chance to have ii(t) = 0. Because the molecules cannot be self-replicated fromnothing the state with n = 0 is an absorbing state: once the system reaches this state there will be no recovery. Usually this will not happen to a species with selective value W > 0, unless it starts with a very small number of molecules (Leung, 1984). In that case the fluctuation will kill it before it has a real chance to grow. It is readily seen from equations (14) and (16) that as t + 00, 4(m) = 1 and u(m) = cc. Therefore the SS will be a stable one with small relative fluctuation d(p/S2). Exactly the same results of (T(M) and P(n; 00) can be derived from the principle of detailed balancing, by having the equilibrium probability function P(n; 00) satisfy the following condition (Van Kampen, 1981), r(n)P(n; 00) = g(n - l)P(n - 1; 00). Furthermore, the previous results (Jones and Leung, well with the present ones if /A < !2 is assumed.

(17) 1981) agree very

4. Stable Selection Resulting from Competition. Within the scheme of the system-size expansion, the average population iii(t) equals the deterministic value G!&(t) for each species in the system. The fate of competition is thus deterministic.

The solution

to equation

(8) is given by,

ABIOMOLECULARSELECTIONMODEL

235

(18) where @i(O) = nj(0)/S2 is the initial concentration. Assuming that W, = max (Wi), we find that as t + OQ,@i(m) = 6ir. Only the superior species with largest W will be selected at the expense of all others. For the stochastic aspect of the selectional processes, we turn to the moment equations for the stochastic components xi(t). Two species with W, > W, will be considered for the sake of simplicity. At the SS predicted for the deterministic components &(t), we have two equations for the first moments, I= .&x;

= 0,

(19)

k

where i = 1, 2. For the second moments, we have the following three,

The SS values of the coefficients that appear in equations (19) and (20) can be evaluated from equations (7) and (12). With t#~:= 1 and 4; = 0, we have,

%!?I&= w, - w,,

$3; = 2A,,

as, = 0.

(21)

From the above, we finally find that 3: = 2; = 0, (J{* = dzz = 0, and

61 =

ccl.

Therefore, a stable selection will be achieved in such a manner that the inferior species are eliminated steadily while the superior one monopolizes the system smoothly with small relative fluctuation ~(~,/s2). The stability can be confirmed by using the linear stability analysis. Comparing with the results of last section, we see that the selection process ends up at the onespecies surviving SS. This agrees with the previous results derived from the moment expansion scheme (Jones and Leung, 198 1). We consider here a special case with all species 5. Metastable Coexistence. having the same selective value W. The deterministic theory predicts a stable SS coexistence with relative concentrations remaining constant in time.

After &(t) achieves the deterministic SS, the behaviors of ail(t) are solved from equation (11) which can be written in a matrix form as,

236

H. K. LEUNG

dS/dt=MS+B.

(23)

For a simple case of two degenerate species with WI = W,= W,S and B are both 3 X 1 column matrices with elements (urr, 022, u12)and (9%&, 0) respectively; M is a 3 X 3 square matrix defined by,

1

(24)

At the SS of Qi, the coefficients in the above are given by, 5Prr =&5

= -w(j”l,

s; = 2/I,@&

&Par=&

= -w@;,

L@;= 2A,9;.

(25)

The SS solution for aft is impossible since det (M) = 0. We thus conclude that after @Qachieves the deterministic SS ori will not be stationary. Their timedependences can be evaluated from equation (23) by using the eigenvalue method. The eigenvalues of M are found to be, Al = 0,

x2=-w,

As=-2W.

(26)

After lengthy but standard. matrix algebra, we find that the magnitudes of oji(t) all contain a common term which increases with time as ~W~~~%P~~~+PZ $S1)at, where

hi = AiIW*

The metastable properties are purely stochastic in nature. The deterministic description of the coexistence is valid only at the early moment of the time-development of ail(t). As will be elaborated in the next section that &(t) = df quickly while oil(t) remains small but increases with time steadily. This metastability can also be described by the probability function P({n}; t) which manages to peak and remain peaked where ni = ii;, and has its width spreading slowly with time. At the later time, the fluctuation will come into play, and the system will eventually become unstable. The possible consequences after the fluctuation catastrophe will be discussed in the next section. We note here that the metabolism factor /.I plays an important role in the stochastic behaviors for the metastable system. Care must be taken in interpreting the results since the validity of the approximation scheme will become questionable when ufj tends to the order of magnitude of #. We cannot here describe the details of the fluctuation catastrophe; instead we predict this tendency. Within the system size expansion scheme, the 6. Stochastic Selection. selectional processes and the final outcomes are purely determinstic:

A BIOMOLECULAR

SELECTION MODEL

231

the species with largest selective value W will be selected to monopolize the system. The stochastic aspects come from the fluctuation effects which depend critically on the metabolism factor cc, and sometimes on the initial populations ~(0). Thus W is not the only parameter to decide the fate of the multispecies competitions. As is discussed in Section 3, ~(0) plays an important role in the early stage of the competition and P is important all the time even after the fittest species is selected since the relative fluctuation of the surviving SS is given byd(p/8). For the special case of having the same W for all species, the deterministic coexistence applies only to {A}, the fluctuation induced transition will occur eventually when P({n}; t) spreads significantly to the absorbing state. We expect that this metastable state will not transit to the absorbing state, instead a stochastic selection will occur in such a way that one species will be selected with the rest eliminated. This conjecture is well justified since we can show easily from equations (8) and ( 11) that the two sums (27) satisfy exactly the stochastic equations for the one species system. The fitness factor will then be p and ~(0) only. The selected species will survive at the one species surviving SS which is both deterministically and stochastitally stable (Leung, 1984). In the case of the completely degenerate system in which all species are identical in terms of W,p and n(O), the selection will be carried out by chance only. The stochastic selection processes are very slow as compared with the deterministic ones which take a time Wt - 10. For the degenerate case, the coexistence in & will also be achieved at Wt - 10, however the lifetime of the metastable state is much longer since Oil(t) are small and increase very slowly. For the completely degenerate case, we may define the lifetime of the metastable state as r = WAt, where At is the time duration starting from the time when {ii} attains the SS values, to the time when the relative fluctuations equal 0.5. We find numerically that, T w O.l116SYJ/~.

(28)

For example, with p = 10 and 52 = 106, we have r - 104. Results derived from the system size expansion scheme agree very well with those from the moment expansion scheme in describing various processes. Both approximation schemes predict the same metastability for the coexistence case; the difference between them is very slight numerically even during the period of increasing ail(t). For instance, the corresponding results of r from the moment expansion scheme is found to be, 7 = 0.1250&?/&

(29)

238

H. K.LEUNG

Both equations (28) and (29) are used here to demonstrate the fact that the lifetime of the metastable state is much larger than the time needed to establish the deterministic coexistence in {ii}. Numerical estimations are surely extended into the region in which the applicabilities of these two schemes are questionable, however the difference between the two is still insignificant. The difference is indeed very small especially at the earlier moment. There is a long time period in which the metastability is described within the accuracy of both approximation schemes. Though we cannot describe the stochastic transition in detail, we can conclude that the predicted metastability will surely lead to the stochastic selection. The author wishes to thank Professor Peter Schuster for his interest and comments regarding the approximation schemes. The partial support of this work by the National Science Council of the Republic of China is also gratefully acknowledged.

LITERATURE of Matter and the Evolution of Biological MacroEigen, M. 1971. “Self-organization molecules.” NaturwiJsenschuften 58,4655 23. and P. Schuster. 1979. The Hypercycle: a Principle of Natural Self-Organization. Berlin: Springer. Gardiner, C. W. 1983. Hundbook of Stochastic Methods. Berlin: Springer. Jones, B. L. and H. K. Leung, 198 1. “Stochastic Analysis of a Nonlinear Model for Selection of Biological Macromolecules.” Bull. math. Biol. 43, 665-680. Kurtz, T. G. 1972. “The Relationship Between Stochastic and Deterministic, Models For Chemical Reactions.” J. Chem. Phys. 57,2976-2978. Leung, H. K. 1984. “Stability Analysis of a Stochastic Model for Biomolecular Selection.” Bull. myth. Biol. 46,399-406. McQuarrie, D. A. 1967. “Stochastic’ Approach to Chemical Kinetics.” J. appl. Rob. 4, 413-478. Oppenhelm, I., K. E. Schuler and G. H. Weiss. 1969. “Stochastic and Deterministic Formulation of Chemical Rate Equations.“J. Chem. Phys. 50, 460466. Van Kampen, N. G. 1976. “The Expansion of the Master Equation.” Adv. Chem. Phys. 34,245309. 198 1. Stochastic Processes in Physics and chemistry. New York. North Holland. -. RECEIVED REVISED

4-13-83 1 l-5-84