On inherited fertility in biological systems: a model of correlated fluctuations in the stochastic branching process

BioSystems, 26 (1992) 185-192

185

Elsevier Scientific Publishers Ireland Ltd.

On inherited fertility in biological systems: a model of correlated fluctuations in the stochastic branching process Yuri A. Berlin, Dmitri O. Drobnitsky, Vitalii I. Goldanskii and Vladimir V. Kuz'min N.N, Semenov Institute of Chemical Physics, Academy of Sciences of USSR, ul. Kossygina 4, 117334 Moscow (Russia) (Received July 2nd, 1991) (Revision received November 22nd, 1991)

A new evolutionary model with hereditary modes considered as correlated fluctuations of fertility has been proposed. It has been demonstrated that the model allows the global statistical properties of the system to be evaluated, e.g. the ensemble average and the probability of extinction. The results obtained show the increase of instability of a population with the enhancement of inheritance efficiency. The existence of at least an exponential stratification in the population has also been shown. Possible applications of the present model are discussed.

Keywords: Correlated fluctuations; Fertility; Heredity; Non-Markovian stochastic branching process.

1. Introduction

There are many biological and physical problems that concern the joint effect of the large number of locally interacting processes invoking each other. These problems appear in the theory of evolution, investigation of populations dynamics, propagation of infections, branching and avalanche physico-chemical processes, etc. (see e.g. Bartlett, 1960; Dremin, 1989; Grossberg and Khohlov, 1987; Koda et al., 1991; Lalley and Sellke, 1987; Ratner et al., 1985; Svirejzev and Logofet, 1978; Tautu, 1988; Watson and Galton, 1874; Zhirmunsky and Kuzmin, 1990). One of the most powerful approaches to the description of such phenomena is based on the model of the stochastic branching process (SBP) (see e.g. Aliev, 1987; Harris, 1963; Karlin, 1968; Moyal, 1962; Watson and Galton, 1874). Various variants of this model have been proposed to describe multi-element biological and physical Correspondence to: Y.A. Berlin.

objects almost since the middle of the last century (Bartlett, 1960; Moyal, 1962; Watson and Galton, 1874) until nowadays (Biggins and G5tz, 1987; Dittrich, 1990, Drobnitsky et al., 1991; Koda et al., 1991, etc.). The great attention paid by investigators to the SBP approach is caused mainly by the dearness of its formulation. Each concrete SBP model is constructed on the local level, i.e. on the level of its elementary act. This allows to take into account characteristics of real systems explicitly. On the other hand, the approach is convenient from the mathematical point of view. It allows conclusions to be made about the global evolution of the system on the basis of assumptions about its local properties. In addition, in most cases the inverse problem can also be solved either directly (see e.g. Zolotarev, 1954) or by means of various fitting procedures (see e.g. Bartlett, 1960; Lalley and Sellke, 1987). The aim of the present work is the further development of the SBP paradigm recalling one of its important features, namely the existence of local fluctuations of reproducing entities

0303-2647/92/$05.00 © 1992 Elsevier Scientific Publishers Ireland Ltd. Printed and Published in Ireland

186

number, i.e. fluctuations taking place in the course of its single act. For instance, such fluctuations can result from variability of individuals in populations. Obviously, in real biological systems, the fluctuations mentioned above may correlate due to environmental and other external causes or due to heredity. This implies that somewhat realistic SBP models should include correlations of fluctuations on the microscopic level. Such a statement of the problem in most cases leads to the violation of the Markovian property in the dynamics of real systems and especially in biological communities (see e.g. Svirejzev and Logofet, 1978; Svirejzev, 1987; Zhirmunsky and Kuzmin, 1990, etc.). Some attempts have been made to consider such correlations in the framework of the SBP model by introducing spatial diffusion of reproducing entities, their mutual transformations, age-dependencies, etc. (Biggins and GStz, 1987; Dittrich, 1990; Harris, 1963; Persson et al., 1991; Rahimov, 1988; Scott and Uhlenbeck, 1942; Yarovaya, 1990 and others). Note that in these cases, local correlated fluctuations of reproducing entities number are not related with innate properties of the entities and virtually constitute only some additional external restrictions on the reproduction process. In contrast to these studies we present here a new variant of the SBP model in which correlated fluctuations exist due to one of such innate properties, namely the inheritance of fertility. Formally our approach constitutes the introduction of hereditary modes into the standard Watson-Galton model (Aliev, 1987; Harris, 1963; Karlin, 1968; Rozanov, 1979; Watson and Galton, 1874). However it can be naturally expanded on the consideration of arbitrary random influences on biological systems. In particular, it allows new insight to be made into such an important problem as the extinction of communities (see e.g. May, 1973; Keiding, 1975; Svirejzev, 1987; Turelli, 1977). The concrete formulation of this approach is given in Section 2. In Section 3, we briefly describe the mathematical aspects of the solution of the problem. The results obtained in the framework of our method are presented and

discussed in Sections 4 and 5. The last section of the paper contains some suggestions about possible generalizations of the method. Finally, we summarize the main conclusions. 2. Model

Consider a number of individuals propagating independently from each other. The number of nearest descendants (children) m of an individual is a random number with a certain distribution ¢(m), m = 1,2,3,.... We specify the mechanism of begetting as follows. All individuals have identical fixed lifetimes ending at the appearance of their children. Thus individuals breed by gemmation. However, the above restrictions on the mechanism of begetting are not crucial for the global behavior of the system, since in the long time limit, the evolution of the system depends only on the average lifetime of individuals and on the distribution ¢(m) (see, e.g. Aliev, 1987; Harris, 1963; Karlin, 1968). Let ~i(t) be the number of descendants (not only nearest ones) of an initially given ith individual which exist at time t. Then the total number of individuals at time t is ~(t) --- ~l(t) + "'" + ~q(O

(1)

given q individuals at t -- 0. If ~(m) is equal for all individuals in all generations, then we are dealing with a simple Watson-Galton branching process which is a Markovian chain in terms of ~(t). We shall call it the process without heredity. Now we introduce heredity in the form gch = f(mp,i')

(2)

where/~ch is the first moment of ~bfor each child of the parent that has yielded mp children. The function f provides the direct correspondence between the expected number of grandchildren and the actual number of children. We shall call f the heredity function. The parameter i" determines the increase or decrease rate of this function and hence designates the form and the

187

strength of correlated fluctuations in the system. Equation (2) implies that children inherit high (or low) fertility of their parents and the productive branches obtain a certain preference if f is a non-decreasing function. The parameter ~ explicates the preference extent. Therefore we shall call the process with the property given by (2) the process with heredity and ~', the factor of preference (FP). For biological and physical systems described by such a process it is important to know the quantities characterizing their main evolutionary trends, namely the extinction probability and the ensemble average of the number of individuals. In principle, this may be done by means of the formalism developed for general branching processes by Harris (1963), Karlin (1968), Moyal (1962) and others. However, its direct application to the present model encounters essential difficulties (Drobnitsky et al., 1991). In the next section we describe the mathematical procedure that calculates the required quantities without referring to the general formalism mentioned above.

and

•

Z ( k , t + 1) =

~k(i) [ZOe(i),t)] i

iffi0

(4) z(k,1)

--

~k(o)

where k is a real positive variable. For the first generation one obviously obtains (X(1)) ffi t~0 (5)

po(1) = ¢~o(0) Recalling the independence of the individuals breeding (see expressions (3)-(5)) and using the induction method, one can directly show that (x(t)) = F(t~o,t)

(6) po(t) -~ Z ( m , t )

given a single initial individual with the feature #0"

3. Mathematical considerations

Recalling the independence of the breeding of individuals, it is sufficient to consider the population with a single initially existing individual with a certain feature (the average productivity) #0. The generalization for an arbitrary number of initial individuals is obvious. Since the process is considered in discrete time, we shall use the lifetime of individuals as the unit of time and mark time in integer numbers. To obtain information about the main evolutionary trends of the system with inherited fertility, we construct the following auxiliary functions: F(k, t + 1) =

~

Ok(i) i F(f(i),t),

iffiO

(8) F(k,1) = k

In contrast to the general formalism (cf. e.g., Harris, 1963), the relations (3) and (4) have appeared to be very convenient for numerical calculations as well as for the analytical evaluations of trend asymptotics and t-dependencies (for more details see Drobnitsky et al. , 1991). As a result the evaluation of the main evolutionary trends of the process with heredity becomes available (see (3)-(6)). 4. Results

To exploit our approach for obtaining numerical results, explicit forms of f and ~b are necessary. These functions are given by the concrete problem. However, in the present paper we are interested only in some typical properties of the system with inherited fertility. Therefore, recalling the illustrative character of the further calculations, we restrict ourselves to the following two, rather representative forms of the heredity function

188

-- rk

(7)

ilk) = kt

(s)

where ~"is a real parameter and k is an integer variable. In addition we assume that ¢ is the Poisson distribution

50-

20"

L1 W

10-

/t m

¢(m) =

m!

exp(- g)

(9)

0

Using the relations (3) and (4), together with (6), for the heredity in the form (7) with ~"> 0, in the long time limit (actually, as valuations show, for t > ~), one obtains (~(t)) - exp(~t 2 + /~t)

B '~

exp('~ q)

(11)

q=O

where ~ - In ~"and/~ - In ~o. For 0 < ~" < 1 as well as for negative L the evaluation of the (~(t)) is available using the expressions (3) and (4), while for ~- > I the relation (11) is valid. Note that the main growth expo-

I

I

,

I

I

I

5

v

,

I

I

~

10

T

I ~-~

15

time(generations)

~ P l ~ s ~

(10)

w h e r e 7 - In~-and/3 - I n g 0 - 1 n L For ~ > 1, the relation (10) provides the growth of the ensemble average in time with the rate greater than exponential. The value ~- -- 1 provides the simple exponential growth (see Fig. la). The main growth exponent ~ increases logarithmically with ~-(see Fig. 2). Note that the simple Watson-Galton process with P[ ~(~) = 01 ¢ 1 provides an exponential growth of the ensemble average with the growth exponent ~ - In t~0 (see Fig. la). For 0 < ~" < 1 the relation (10) provides the decrease of the ensemble average down to zero for t - ~ , while for short times an increase of the ensemble average may occur (see Fig. lc). In this case the increase and decrease rates are greater than exponential due to (10). Using the relations (3) and (4) together with (6) for the heredity in the form (8) with ~- > 1, one obtains

(}(t)} - exp

1

0

~ -'~

L~(~.~)

P2

time c

.._

L~

tt~

:,-,,, 0

\ 5

\ 10

time (generations)

15

Fig. 1. The growth of the ensemble average of the number of individuals. (a) L1, heredity in the form (7) with ~- = 1; L2, heredity in the form (7) with ]" > 1; P, heredity in the form (8) with ~ > 1; W, Watson-Galton process. All dependencies are plotted in the logarithmic scale. (b) P1, heredity in the form (8) with 0 < ~" < 1; the dashed line is the asymptote; P2, heredity in the form (8) with ~" = 0. The dependencies are plotted in the logarithmic scale. (c) L1 and L2, heredity in the form (7) with 0 < ~" < 1 (~L1 > ~'L2);P, heredity in the form (8) with ~" < 0.

189

9,

3.0 -~

po(O

p: 2.0

1.0 .... i---'

0 0.0

0.0

..... 2.0

7*

I

4.0

6.0

i

I

'

~-

--r

....

r - ~ 7

......

r----T"--7

....

5 bim e (g ~ n e r a Novrs )

r-

-7

10

I

8.0 10.0

¢

Fig. 3. The growth of the extinction probability with time (qualitative picture).

Fig. 2. The dependence of the main growth exponent 7 versus the FP value (~).

nent also increases logarithmically with i" as in the previous case (Fig. 2). Thus for ~ > 1 the ensemble average of the process with heredity in the form (8) demonstrates an overdrive (exponent versus exponent) growth (see Fig. la). For ~" = 1, equation (8) reduces to the case considered above (ilk) = k). For 0 < ~ < 1 the calculation due to (3) and (4) shows that (~(t)) grows slower than an exponent but reaches the exponential growth asymptotically (see Fig. lb). For ~" = 0, one arrives at the set of trivial W a t s o n - G a l t o n processes (with ~ constant 1) starting from the second generation. In this case (~(t)) = constant = g0 (see Fig. lb). For ~ < 0, the heredity functionJ(k,D with k = 0 should be redefined. Defining j ( 0 , ~ ) = 0, one obtains the overdrive decrease of (~(t)) (see Fig. lb). The extinction probability, po(t) =- P{~(t) = 0 }, obviously increases with time and reaches its asymptote Po - p{~(co) = 0} as shown in Fig. 3. Since the probability distribution p,(t) =P{~(t) = n} is defined on positive integer numbers n, it can be shown that po(t) = 1 when and only when (~(t)) = 0. Hence in all cases considered above, providing (~(t - co)) = 0, one obtains P0 = 1. In other cases 0 < p0 < i if ~ ( 0 ) > O. Using the relations (4), it can be directly shown that p 0 decreases with the increase of (if the increase of ~"leads to the reinforcement of the heredity function growth rate as it is in our =

=

examples). Furthermore, i f f ( k , t - Qo) - co for any k > 0, then P0 - (¢~0(0) and i f f ( k , t - oo) constant < Qo for some k > 0, then Po constant > ~b~o(0). For instance, for the heredity in the form (8), one obtains 0 < ¢~o(0) < p o < 1 (see Fig. 4). For the heredity in the form (7), Po converges to ¢~0(0) with increase in i" (see Fig. 4). An interesting feature of the heredity function (7) with 0 < ~ < 1 is the relay of the increase and the decrease of (~(t)). Such peculiar behavior of the system is determined by the production of the large number of low-productive individuals in first generations. Initially they cause the maximum of (~(t)) at some time t, but they yield few descendants which, in addition, have low pro-

po*

a>#MO) a=#~Co)

( Fig. 4. Reaching the asymptotes (dashed line) a by the full extinction probabilities Po with increasing j'-value. L, heredity in the form (7); P, heredity in the form (8). Both are qualitative pictures.

190

ductivity. That is why in the following, the population becomes extinct. Summarizing the results, we can conclude that the existence of inherited fertility makes the evolution of the system more crucial, i.e. the greater is FP (represented by the i" value), the less time is required to determine the fate of the system with respect to growth-extinction trends. 5. Discussion

5.1. Non-Markovity in evolution It is becoming increasingly evident that the consideration of random and regular impacts on populations and biological communities leads to the essentially non-Markovian character of populational dynamics as well as biological evolution in general. Therefore statements and solutions of non-Markovian evolutionary problems are of great significance for modern biology and ecology. In the present work, the violation of the Markovian property is caused by correlated fluctuations appearing due to the inheritance of fertility. As follows from the results, evolution of the population with the inherited fertility demonstrates much more interesting behavior than the simple Watson-Galton model (cf. e.g. Aliev, 1987; Bartlett, 1960; Harris, 1963; Watson and Galton, 1874). In particular, increase in the population instability with the reinforcement of inheritance efficiency has been shown.

certain probability. The latter changes with time due to future mixing. Every individual produces a variety of genetic signs in each begetting act. Hence the global variety is proportional to the ensemble average. Thus mixing takes place together with the stratification of individuals due to their fertility. As is clear from the results, while the average increases at least exponentially in time, the extinction probability remains approximately constant. The set of probabilities {P{~(t) = n l}~ = 0 is obviously normalized at any fixed time instant t. Consequently at large times this set constitutes a double-peak distribution as shown qualitatively in Fig. 5. The same situation appears in processes without heredity. However in the process with heredity the right-hand peak runs forward at least exponentially faster (Fig. 5). Thus mixing together with stratification take place in the exponential regime on average.

5.3. Inverse problem An important consequence of our consideration is the possibility of evaluating the information concerning the hereditary mode from the experimental data. Such a possibility appears if an observer is able to change the F P value. Then using ~(~) and p*(~') as well as (~(t)) and po(t) discussed in Section 4, one can make a conclusion about the heredity function by fitting. The latter contains all the individual's genetic information.

5.2. Mixing and stratification An interesting feature of the model is the permanent perplexing of genetic signs, i.e the system starting from the pure set {~i} evolves to a mixed state. The sparse genetic signs (parameters of ~) caused by correlated fluctuations increase together with the ensemble average. This leads to the principal impossibility for an observer to extract pure genetic signs from such a final state. The only conclusion we can make about the particular genetic sign in principle is that it exists in the mixture with a

f.

'

/4

Fig. 5. The qualitative form of the probability distribution p~(t) ~ P I $(t) = n} in the long time limit for SBPs, provid-

ing the growth of the ensemble average; solid line, the process with heredity (~ > 1); dashed line, the Watson-Galton process. The arrows mark the direction of the movement of the right-hand peaks.

191

5.4. Intensity and: productivity

tion probability of the group of cells with fluctuating productivity will be comparatively large.

In our model the intensity of gemmation (the lifetime of individuals) is assumed to be fixed, and the fluctuating value is fertility (set {~b(m)}). But the transition to the heredity in terms of intensity is straightforward. Assume that the set ~b(m)is fixed and that the lifetimes of the individuals are random and are distributed due to the exponential expectation law. Let r -- 1/k be the average lifetime of an individual. Then (2) can be replaced by kch -- g(kp)

(12)

Defining local transition densities (see e.g. Karlin, 1968; RozsLnov, 1979) as k~ = k~b(m)

(13)

we obtain the fluctuating set (kin) defining the reproduction. However the same ks can be obtained by assuming k = constant equal to the initial value and tLch = flmp) and by matching f and g using (12) and (13), i.e. deforming ¢(k) instead of varying ~. Note that in this case we must measure time in units of initial 1/k since the latter is now the unique scale factor. Deforming and matching procedure allow the problem to be solved using the formalism considered in Section 3. As a result the proposed approach can be applied not only to the processes with inherited fertility but also to those with inherited intensity. For example, let us consider the cell division with ¢(m) ~ 0 for m = 0,1,2 and Pm = 0 for m > 2. In addition, suppose that the inheriting feature is the division intensity. Replacing g b y f w e obtain some new set {~(m)}. Obviously, ~b(m) would not necessarily be the Poisson distribution. However one can expect that the behavior of the population will be similar to that described in Section 4. The system will reach its asymptotics (will produce a large number of cells due to productive ones or will collapse) after 4 - 8 normal generations. The intensively dividing cells will have yielded probably hundreds of generations by this time. However the extinc-

6. Conclusion

The model presented concerns the mechanism of the selection of states due to a certain feature. We have found the main features of the evolutionary behavior of systems with the feature controlled by the factor of preference. This preference is created by correlated fluctuations postulated in the model and described in terms of the heredity function. Such a function can be evaluated from experimental data by fitting procedures. For more complicated forms of heredity, equations (2) and (12) should be replaced by general functional equations. However if the fertility is inherited, the main results are likely to be similar to those obtained in Section 4. Thus the appearance of inherited productivity in the form presented in biological systems leads to a certain sharpening of its two attractors (disappearance and explosion, see Fig. 5) and hence to the increase of instability. Situations similar to those considered in Sections 4 and 5, and in particular those with nontrivial behavior, appear in various population models taking into account the influence of the provision of resources, concurrent struggle, environmental conditions, etc. (see e.g. Bartlett, 1960; Coste et al., 1978; Dawson, 1975; Dittrich, 1990; Eigen et al., 1988; Keiding, 1977; May, 1973; Moyal, 1962; Ratner et al., 1985; Svirejzev and Logofet, 1978; Tautu, 1988; Turelli, 1977; Yarovaya, 1990; Zhirmunsky and Kuzmin, 1990, etc.). This implies that the explicit form of correlated fluctuations in the system given by the relation (2) or (13) can be interpreted not only in terms of intrinsic properties of individuals but also in terms of the joint effect of innate and environmental impacts. This fact constitutes the evidence in favor of the general nature of the model presented above. The model considered is applicable not only to population dynamics. Various branching biological, physico-chemical and multi-particle production processes may include correlations

192

of local productivity fluctuations. Our approach and its possible generalizations can provide useful methods for investigation of such phenomena. References Aliev, S. and Shirenkov, V., 1987, Asymptotic Problems of Stochastic Processes Theory. An Asymptotic Behavior of Watson-Galton Processes Close to Critical (Nauka, Kiev), pp. 5,6 (in Russian). Barlett, M.S., 1960, Stochastic Population Models in Ecology and Epidemiology (Methuen, London) pp. 10- 72. Biggins, J.D. and GStz, Th., 1987, Expected population size in the generation-dependent branching process. J. Appl. Prob. 24, 304-314. Coste, J. et al., 1978, About the theory of competing species. Theor. Pop. Biol. 14, 165-184. Dittrich, P., 1990, A critical branching process in random environment. Theory Prob. Appl. 35, 612-615. Dawson, D.A., 1975, Stochastic evolution equations and related measures processes. J. Multivar. Anal. 5, 505 - 508. Drobnitsky et al., 1991, Correlated fluctuations in multielement systems: the stochastic branching process model. Phys. Rev. A, in press. Dremin, I.M., 1989, Correlations and fluctuations in multipartical production. Usp. Fiz. Nauk, 160, 105 - 133 (in Russian). Eigen, M. et al., 1988 Molecular quasi-species. J. Phys. Chem. 92, 6881-6891. Grossberg, A.Yu. and Khohlov, A.R., 1987, Statistical Physics of Macro-Molecules (Nauka, Moscow) (in Russian). Harris, T.E., 1963, The theory of branching processes, in: Die Grundlagen der Mathematischen Wissenschaften band 119 (Springer-Verlag, Berlin)pp. 11-350. Karlin, S.A., 1968, A First Course in Stochastic Processes (Academic Press, New York). Keiding, N., 1975, Extinction and exponential growth in a random environment. Theor. Pop. Biol. 5, 117-158. Lalley, S. and Sellke, T. 1987, A conditional limit theorem for the frontier of branching Brownian motion. Ann. Prob. 15, 1052-1061.

Koda, S. et al., 1991, Branching ratios in O(P - 3) reactions of terminal olefins studied by kinetic microwave absorption spectroscopy. J. Phys. Chem. 95, 1241-1244. May, R.M., 1973, Stability and Complexity in Model Ecosystems (Princeton University Press, Princeton, N J). Moyal, J.E., 1962, Multiplicative population chains. Proc. R. Soc. London Ser. A 226, 518-526. Nummelin, E., 1984, General Irreducible Markov Chains and Non-Negative Operators (Cambridge University Press}. Person, M.D. et al., 1991, The influence of parent bending motion on branching at a conical intersection in the photodissociation of CH31, CD31, CF31. J. Chem. Phys. 94, 2557- 2563. Rahimov, I., 1988, Local limit theorems for the critical Watson-Galton processes with deceasing immigration. Theory Prob. Appl. 33, 387-392. Ratner, V.A. et al., 1985, Problems of Molecular Evolution Theory (Nauka, Novosibirsk} (in Russian}. Rozanov, Y., 1979, Random Processes (Nauka, Moscow} (in Russian}. Scott, W.T. and Uhlenbeck, G.E, 1942, On the theory of cosmic-ray showers. Phys. Rev. 62, 497-508. Svirejzev, Yu. and Logofet, D., 1978, On Stability of Biological Objects (Nauka, Moscow) pp. 119-123 (in Russian). Svirejzev, Yu., 1987, Nonlinear Waves, Dissipative Structures and Catastrophes in Ecology (Nauka, Moscow) pp. 352- 355 (in Russian). Tautu, P., 1988, On the qualitative behavior of interacting biological cell systems. Stochastic Processes in Physics and Engineering, Dordrecht, pp. 381-402. Turelli, M., 1977, Random environment and stochastic calculus. Theor. Pop. Biol. 7, 140-172. Watson, H.W. and Galton, F. 1874, On the probability of extinction of families. J. Anthropol. Inst. Great Britain and Ireland. 4, 138-144. Yarovaya, E., 1990, Alteration of high order moments in the model of branching process with the diffusion in a random media. Vestnik MGU 1, 79-82 (in Russian). Zhirmunsky, A.V. and Kuzmin, V.I., 1990, Critical Levels in the Developement of Natural Systems (Nauka, Leningrad) (in Russian). Zolotaryev, V., 1954, One inverse problem in the theory of stochastic branching processes. Usp. Mat. Nauk. 9, 147-156 (in Russian).