key: cord-0832410-qmp8g39u
authors: Swartz, Daniel; Ottino-Loffler, Bertrand; Kardar, Mehran
title: A Seascape Origin of Richards Growth
date: 2021-08-23
journal: Physical review. E
DOI: 10.1103/physreve.105.014417
sha: c7ed6e6d29c4696c8f8dc55ee59f0103ce879c19
doc_id: 832410
cord_uid: qmp8g39u

First proposed as an empirical rule over half a century ago, the Richards growth equation has been frequently invoked in population modeling and pandemic forecasting. Central to this model is the advent of a fractional exponent $gamma$, typically fitted to the data. While various motivations for this non-analytical form have been proposed, it is still considered foremost an empirical fitting procedure. Here, we find that Richards-like growth laws emerge naturally from generic analytical growth rules in a distributed population, upon inclusion of {bf (i)} migration (spatial diffusion) amongst different locales, and {bf (ii)} stochasticity in the growth rate, also known as"seascape noise."The latter leads to a wide (power-law) distribution in local population number that, while smoothened through the former, can still result in a fractional growth law for the overall population. This justification of the Richards growth law thus provides a testable connection to the distribution of constituents of the population.

The mathematical description of growing populations has been an enduring topic of interest. There is presently a great diversity of population growth models [1] [2] [3] , which have been applied in a wide variety of contexts, including epidemiology [4, 5] , forestry [6] , developmental biology [7] , cancer [2] , immunology [8] [9] [10] [11] , and many others.

One such model of interest is the Richards Equation. Originally developed as an empirical model of plant development [6, [12] [13] [14] , and later used as a population growth model [15] , it has over the last year exploded in popularity, becoming a consistent presence in the modeling of infectious diseases, especially in COVID-19 forecasting [4, [16] [17] [18] [19] [20] [21] [22] [23] . The Richards model is described by

where y is a measure of the population size, µ is the initial growth rate (sometimes referred to as the population fitness), and a sets the final (saturation) population to y f = (µ/a) 1 γ−1 . The feature that distinguishes this from the celebrated and more natural Verhulst logistic equation [24] is the non-analytic (shape) parameter γ > 1. When γ = 2, a logistic equation is retrieved, but for any other value, the inflection point of the growth curve becomes off-center, leading to asymmetric growth curves, as shown in Fig. 1 .

This shape parameter γ indeed gives the Richards equation character, but its meaning beyond a fitting device is mysterious. Because γ takes on factional values, this makes the dynamics of Richards growth nonanalytic, making its origins theoretically nontrivial. Although there have been proposed derivations of Richards-like growth, they have either leaned on an underlying fractal topology [19, 25, 26] , or relied on a detailed manipulation of an SIR model, which may not be robust under model perturbation [18, 27] . In practice, Richards growth is still considered an empirical law and lacks a fully universal origin, or a physically intuitive interpretation of the shape parameter γ.

In this paper, we propose a plausible origin for emergence of Richards-like growth in distributed populations from generic local analytical forms. Specifically, we shall use the Fisher Equation as starting point. The deterministic Fisher equation [28] [29] [30] [31] [32] is one of the most basic models of spatial population growth, written as dy(x, t) dt = µy − ay 2 + D∇ 2 y .

This is distinguished from a logistic equation by the presence of a spatial coordinate x and a diffusion term D which sets the rate of local migration. The deterministic Eq. (2) should in principle include reproductive stochasticity. While there are many kinds of randomness that may influence population growth (e.g., demographic noise, with amplitude proportional to √ y [33] [34] [35] [36] [37] [38] ), we are interested in what is called Seascape Noise [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] . This noise arises from observing that the fitness µ of a population is not necessarily a static, uniform quantity. Population fitness can vary based on minor environmental conditions, access to resources, microscopic mutation rates, and other such factors. So we generalize from a static fitness landscape µ(x) to one that varies in time µ(x, t), much as the sea surface changes in time. By introducing seascape noise into the Fisher equation, we obtain the Seascape Fisher Equation

with σ 2 being the variance of the fitness noise, and dW = dW (x, t) is an uncorrelated Wiener process. While our noise will be uncorrelated, we should note certain applications may wish to consider more spatially correlated environmental variations [49] [50] [51] [52] . Moreover, we consider this stochastic differential equation under an Itô interpretation, since its assumptions pair well with the fact that generations in a popuation are discrete, though interesting seascape results exist for alternate interpretations [53] [54] [55] .

Since both the Fisher equation and seascape noise are very versatile concepts, this makes their combination in Eq. (3) similarly adaptable. Variants of this much-studied model [56] [57] [58] [59] have connections to the Kardar-Parisi-Zhang equation in the field of surface growth problems [60] [61] [62] [63] [64] [65] [66] [67] [68] , the study of directed polymers in random media [56] [57] [58] [59] , symmetry breaking [69] [70] [71] , theoretical ecology [35, [72] [73] [74] [75] [76] , immunology [77] , and economics [78] .

In this work, we argue that spatially averaging the seascape Fisher equation is sufficient to produce Richards-like growth curves. In Sec. II we introduce a discretization of the model, replacing y(x) with y i for locations (nodes) i = 1, . . . , N . In the context of this model diffusion between neighboring sites can be replaced with migration between any pair of sites at rates M ij . A mean-field limit is then obtained if migration rates are equal between any pair of sites (M ij = D/N ). This limit was considered by us [79] in the context of population extinction, and similar models by others in different contexts [69, 70, 77, 80, 81] . The steady steady (long-time) limit of the probability distribution ρ(y) can be obtained exactly in this meanfield limit, and characterized by a power-law distribution. A wide power-law distribution, characterized by non-analytic dependencies of its moments, suggests a natural route for obtaining Richards-like growth. Unfortunately, we could not solve the full dynamic behavior of the model, even in the meanfield limit. As an alternative, in Sec. III we introduce a model which alternates the linear and nonlinear parts of the model. This seasonal growth model retains the features of seascape stochasticity and migration which we believe are the cause of Richards growth, and (with some assumptions) is exactly solvable; numerical simulations confirm these expectations. In Sec. IV we consider a number of extensions: providing a procedure for testing the model by examining the dependence of the second moment on the mean (Sec. IV A); indicating the universality of the results for generalized growth equations (Sec. IV B); numerical studies of the seascape Fisher equation in one and two dimensions (Sec. IV C); and finally indicating the appearance of Gompertz growth in the strong noise limit (Sec. IV D). The concluding Sec. V provides an overview, and outlook for future studies.

The intuition for why seascape noise is a viable mechanism for Richards growth comes from earlier studies of a mean-field version of the problem. We begin with the full seascape Fisher equation Eq. (3) in d spatial dimensions. We discretize space into N lattice sites with unit spacing indexed by i = 1, 2, . . . , N − 1 with periodic boundary conditions. The field y(x, t) is replaced by a single value at each lattice site, y(x = i, t) → y i (t). The Laplacian on this lattice takes the form of a migration rate M ij giving the rate of migration from site j to i. For example, on a regular one-dimensional lattice, M ij = δ i+1,j + δ i−1,j − 2δ ij . The discretized version of the stochastic Fisher equation now takes the form

One can use the above equation and generalize the matrix M ij to any migration connectivity network. We then proceed to take the mean field limit by using a complete graph, defining M ij = D/N , cor-responding to a population that can migrate with "steps" of arbitrary length. The stochastic Fisher equation in the mean-field case now takes the simple form

whereȳ is the spatial average of y [82] . The steady state (long-time) behavior of this stochastic ordinary differential equation is obtained as the stationary solution of a corresponding Fokker-Planck equation (see Ref. [79] and Sec. III A), and is proportional to

where c D = 2D/σ 2 , c µ = 2µ/σ 2 , c a = 2a/σ 2 . In the limit N → ∞ we can identifyȳ = y as a parameter to be found self-consistently [69, 70, [79] [80] [81] 83] . The steady state distribution thus has a power-law with upper and lower cutoffs of c −1 a and c Dȳ respectively. Distributions of this form give rise to anomalous scaling in the moments, which we conjecture as the mechanism for Richards growth.

Demanding thatȳ be the mean of this distribution implies the self consistency condition y = yρ(y) dy ρ(y) dy .

By manipulating this condition [79] , we arrive at a moment scaling relationship of

Thus, for c D − c µ < 1 we have that the second moment scales with the first moment as y 2 ∼ y 1+c D −cµ , while for c D − c µ > 1, we have a more typical scaling y 2 ∼ȳ 2 . This anomalous scaling in steady state is what motivated our hypothesis for Richards growth. However, absent a solution of the time-dependent Fokker-Planck equation, we cannot identify parameter regimes where the previous scaling holds away from steady state. As such, we introduce a seasonal model that allows us to take advantage of the results of this section. A justification for employing the steady-state solution in a dynamic setting comes from the case of population decay for µ = 0: If we take the scaling, y 2 ∝ y 1+min(c D ,1) , from the stationary-state, and assume that it is maintained quasi-statically during the decay process, we can make a prediction for the decay of the mean. In particular, ∂ t y = −a y 2 ∝ − y 1+c D , indicates that the mean should scale as t −1/c D . The numerical results plotted in Fig. 2 , indicate a power-law decay roughly consistent with the this quasi-static approximation. A similar argument is made in Ref. [84] where power-law decays as in the quasi-static assumption during growth, we next introduce a two-state model next that allows to more controllably take advantage of this solution.

The combination of nonlinearity and stochasticity complicates the study of Eqs. (3) and (5). To simplify analysis while maintaining the qualitative features we believe to be responsible for the Richards form, we separate non-linear growth and stochastic migration. To do so, we introduce a seasonal growth model composed of two distinct stages:

1. Exploration: The exploration phase (season) has the population diffusing in the presence of seascape noise, as described by

In the mean-field limit, D∇ 2 y is replaced by D(ȳ − y), whereȳ = y at the start of the exploration phase. The exploration phase is run for a time T e with DT e 1.

As these dynamics conserve the mean,ȳ can be treated as a constant. In the absence of reproduction in the exploratory phase, a new interpretation for seascape noise is necessary. A potential cause is random extinctions of the populations at different locales and times. In this interpretation, the noise must have a negative mean, and y will decrease through the process. For a large enough population (to avoid the possibility of extinction) this overall loss can be restored in the subsequent growth stage, without qualitatively changing the results.

2. Growth: Even when starting with a uniform distribution of numbers at each node, the stochasticity in the exploratory phase leads to a broad distribution of y i at the end of this interval. In the subsequent growth phase, at each node we implement reproductive growth following the logistic equation,

Equation (2) is deterministic, and can easily be solved to give the final population at each node in terms of the initial value. The exact form of the equation governing growth and saturation is not important, and as discussed in the following, any analytical form leads to similar results.

We alternate between exploration and growth for a large number of times (seasons) to generate trajectories for the moments of the population m n (t) = y(t) n .

In the exploration phase of the mean-field model, the distribution of y evolves according to the Fokker-Planck equation

The stationary steady state solution is found by demanding a vanishing probability flux, and is proportional to

which has the advertised structure of a power-law with c D = 2D/σ 2 , and a lower cutoff set by c Dȳ . Power laws arising from models with seascape noise are quite common, so this is well-precedented [8, 35, 77, [85] [86] [87] [88] . However, the continuous variation of the exponent with the ratio c D is unusual and likely a feature of the mean-field limit; though some similarly continuous non-universal exponents have also appeared in directed percolation under temporal disorder [89, 90] . However, the steady state solution is only valid in the limit T e → ∞; for any finite T e , the power-law tail won't extend to infinity, and higher moments will not diverge (see Appendix C). In addition, expansion and growth occur simultaneously in the original seascape Fisher equation, making the establishment of the full power-law tail implausible. Because of this, we also impose an upper cutoff Λ to this distribution, simplifying it to

with ρ(y) = 0 otherwise. We have additionally imposed a lower cutoff Υ to ensure that the mean is equal toȳ. For 0 < c D < 1, and as long as Λ Υ, the first moment of the distribution is given by

setting the lower cutoff to

The second moment is now given by

B. Deterministic growth

Once the growth phase begins, the variable y i of nodes, distributed according to Eq. (11), serve as initial conditions for the logistic equation. For each node, the solution of the logistic equation gives y i (t) = e µt µy i ay i (e µt − 1) + µ .

Averaging over the initial conditions now leads to

While this integral can be approximated in detail (see Appendix B), as long as the growth phase only occurs for a small time (µT g 1) we can expand our result for the mean at the end of the growth phase to O(T g ) and find y(t + ∆t) = (1 + µ∆t) y(t) − a y 2 (t) ∆t . (17) If we take advantage of the fact that Λ Υ, then we can take advantage of the previously computed moments to obtain . Therefore, when we write the dynamics of the mean in the growth phase, we have one term that scales as y and one term that scales as y 1+c D , leading to non-analytic, Richards-like growth.

Upon iterating these dynamics at times t j = j(T e + T g ), we find a population trajectory m 1 (t j ) = y(t j ) governed by approximate Richards law to leading order in time andȳ/Λ, of

whereμ = 1 1+Te/Tg µ is the effective growth rate and is cutoff-independent. The carrying capacity is set byã(Λ) which is some function of the cutoff and can be found numerically. We confirm these results in Fig. 3 where we see excellent agreement between the simulated mean population and the analytical prediction. We also note that the numerical results can also be reasonably fitted to a logistic sigmoid. However, what is significant is that models such as Eq. (3) contain the essential ingredients to analytically justify a Richards growth law.

In this section we test the generality and applicability of the results from the seasonal model in a number of other contexts.

In the seasonal model the non-analytic form of Richards growth emerged from the power-law character of the distribution of local numbers, established from the combination of migration and seascape noise in the exploration phase. In the original model, and many likely realizations, growth and migration occur simultaneously, preventing formation of a well defined quasi-stationary power-law distribution. This complicated testing the foundation of our explanation of Richard's growth in terms of a relation between the growth exponent of the distributed population and its distribution amongst different locations.

To address these questions we simulated the full dynamics of Eq. (5) for N = 2 20 nodes. The mean, y = m 1 , of this distributed population grows in sigmoid fashion as depicted in Fig. 4a . In the absence of noise, with µ = a = 1, the population would have saturated to µ/a = 1; seascape noise reduces the saturation value to roughly 0.64. Allowing for a variable Richards growth exponent, the sigmoid curve is well-fitted with a value of γ = 1.79. (Note that with the simultaneous action of growth and diffusion, we do not expect γ = 1 + c D , with c D = 1.2 for the chosen parameters of σ 2 = 2 and D = 1.2.)

Nonetheless, we can test whether something like our proposed mechanism is at work by plotting the time evolution of the second second moment y 2 = m 2 as a function of the first y = m 1 . As depicted in Fig. 4b , at intermediate and late times, we find m 2 ∝ m 2 γ -a straight line in the log-log plot with a slope of γ = 1.79. Such log-log fits are a convenient way to retrieve the growth exponent γ in models such as these, and will be used throughout this paper. (The intercept for the log-log line is also useful, and related to the carrying capacity.)

Finding the first and second moments of a distributed population in a natural setting should be easier than characterizing the whole distribution. A relation such as m 2 ∝ m 2 γ , coupled with a Richards fit to the growth curve with the same exponent would provide a good test of the proposed hypothesis.

The formulation of the model in terms of the logistic equation may give the impression that the obtained results are a consequence of its quadratic form. There are in fact many other analytic growth functions with saturation that generate sigmoid curves. The anomalous scaling, induced by the heavy-tailed distributions attributed to spatial dif- fusion and seascape noise, actually confers a universality to these results, independent of the assumed growth equation. Let us consider a generic analytical growth curve g(y) with an initial growth rate µ that saturates to a carrying capacity y = K, which can be written as

The analytic function f (y), with f (0) = 1, alters the growth away from the purely logistic case, and we will demand that f (y) > 0 for all y to ensure the only fixed points of the local growth are at y = 0 and y = K. This also guarantees that our local dynamics given g(y) only have simple (non-repeated) roots. We will again demonstrate how such a generalized model behaves in the seasonal dynamics introduced before, where the population alternates dispersal with seascape noise, and local growth according to Eq. (20) . As before, we restrict analysis to mean-field dispersion during the exploration phase, leading to

Growth . numbers, as in Eq. (11) . With the mean set by y = y ∝ Υ, (for Λ Υ) the higher moments behave as

(22) Notably, each higher moment has the same scaling withȳ, although the prefactors vary by powers of the upper cutoff Λ.

Following a localized growth step of duration ∆t µ −1 , the mean population is increased by ∆t g(y) . Since g(y) can be expanded as a power series g(y) = µy + g 2 y 2 + g 3 y 3 + · · · , we obtain

(23) Once again, averaging over the exploration distribution establishes a Richards growth law; the difference with the purely logistic case is absorbed into the (cut-off dependent) coefficient of the non-analytic saturation term. We test this universality by simulating seasonal growth with a quadratic f (y) curve in Fig. 5 .

Armed with insights from the mean-field model, we now return to the seascape Fisher Eq. so, we numerically simulate Eq. (3) using a centraldifference discretization of the Laplacian on a lattice. We observe that the mean population can be well fit by a Richards curve, with the Richards exponent and carrying capacity are numerically determined from the data as in the mean-field case (Sec. IV A). As depicted in Fig. 6 this procedure leads to a characterization of the one dimensional growth with an emergent Richards shape parameter of γ ≈ 1.75. The corresponding results in two dimensions (for a square lattice and periodic boundary conditions) are presented in Fig. 7 with shape parameter γ ≈ 1.71. Here, we present these shape parameters as effective exponents describing particular simulations. Further investigation, possibly following Refs. [58, 84, 91] may shed light on the values of these exponent and their universality.

We can rewrite Eq. (18) (see also Eq. (23)) as

in terms of a final capacity K. In the limit of strong seascape noise, we have c D = 2D/σ 2 → 0, and if we 

The above form, known as the Gompertz equation is another commonly used growth law [1, 2, 13] . The qualitative prediction of the model is that the Gompertz law should emerge in systems subject to strong seascape noise.

Growth and saturation phenomena are prevalent in nature, leads a host of models and mathematical forms to quantify their description. The Richards and Gompertz laws stand prominently by providing successful empirical fits to sigmoid curves in diverse contexts [1, 2, 4, 6, 7, 12-15, 19-21, 26] . A typical application involves data obtained by summing (averaging) numbers from a distributed dataset. The appearance and success of such non-analytical forms in describing sums is quite surprising: In the same sense that the probability distribution for the sum of many variables typically assumes the analytic Gaussian form, the most natural time evolution for the sum is an analytic growth law (whose first two terms form the logistic equation). A non-analytic series (much like a critical exponent in critical phenomena) requires specific justification.

In this work, we argue that there is good theoretical and numerical evidence to suggest that Richardslike growth emerges as a natural consequence of the combination of dispersion and (seascape) stochasticity in a large distributed population. In steady state, the combination of the two establishes a broad power-law distribution at distinct locales; averaging any analytic growth rule with such a distribution leads to a non-analytic Richards form. This argument, which nicely connects distribution of local numbers to the time evolution of the mean, however relies on a form of quasistatic evolution that cannot be rigorously justified except in the case of a seasonal model that separates growth and stochastic dispersion. Fortunately, our numerical results in a number cases other than the seasonal model suggest the broader applicability of this result. We further suggest a scheme (plotting the second moment as a function of the mean) to test the validity of our proposed mechanism.

Several aspects of the model that require further study: While the quasistatic arguments lead to a shape exponent of γ = 1 + 2D/σ 2 , this is not the observed exponent even in the mean-field limit with a finite µ. The complexities of an evolving probability distribution invalidate the quasistatic result. It may be possible that a more delicate treatment of the Fokker-Planck equation can help identify an effective Richards exponent. The study of growth laws, or even statics, for the seascape Fisher equa-tion in finite dimensions is also interesting, and may have connections to directed percolation [92] and directed polymers in random media [93, 94] .

More complex connectivities are appropriate to social networks and spreading of contagions. Considering how the Richards equation has seen a recent rebirth in epidemiology and pandemic forecasting, it would be interesting to see whether implementing seascape growth on such human networks would generate relevant results. In particular, we may explore the role of long-range interactions may be important to Richards growth [25, 26] .

While an experimental model can have complicating features, a bacterial population in a well mixed fluid may present a controllable realization of the mean-field model. Perhaps local growth rates can be varied stochastically by application of randomized light sources. The prediction of the model is that upon increasing the strength of reproductive noise, the overall growth curve will change from logistic to Richards to Gompertz. which, if we neglect the higher order terms, returns us to Eq. (18) . Therefore, a Richards-like growth law is retrieved.

If we are interested in the limit of small c D , notice how the first term of the expansion has a term like (1 − x c D )/c D . Therefore, the limit becomes

which is difficult to rewrite in terms of y alone.

Appendix C: Moments of Seascape Dispersion while not directly relevant to the study of scaling laws, for completeness we include a study of the dynamics of the moments, which can be used to determine relevant timescales for simulation. Let's consider a system under seascape noise and dispersion, much like the exploratory phase of the seasonal model of Sec. III, with dy = D(ȳ − y)dt + σydW .

(C1)

By applying Itô's lemma and taking the expectation of both sides, we get a set of ordinary differential equations of the form Note that each moment m depends only on itself and the moment m − 1, meaning that this is a solvable system. For example, because y(t) =ȳ for all time, this implies that the second moment is given by

where h m ≡ y m (0) , and assuming generic values of σ and D (such that the denominator is never 0). Notably, for c D = 2D/σ 2 < 1, this means that the second moment grows exponentially. Given that the (m − 1)'th moment is known, the m'th moment is given by 

The coefficients C m,k can be solved for exactly using inductive methods, though such precision is mostly of specialized use. More importantly, this relation means that each moment above m = 1 experiences leading-order exponential growth, with a growth rate of m[(m − 1)σ 2 /2 − D] > 0 since c D < 1. Using this information, we can reasonably choose a befitting timescale for the exploration phase.

Infectious Disease Modelling

Theoretical population biology

The adaptive seascape: the mechanism of evolution

Advances in Condensed Matter and Statistical Physics

Physical review letters

Advances in physics

Proceedings of the National Academy of

Journal of nonlinear science

Becauseȳ depends on the solution ρ, this makes the SDE (5) a McKean-Vlasov equation

Proceedings of the National Academy of Sciences of the United States of America

Ecole d'été de probabilités de Saint-Flour XIX-1989

Stochastic processes and their applications

We describe here the numerical tools used to obtain the simulation results. The integration of stochastic equations is in general more difficult than the deterministic case. Deterministic ODEs can be treated with a variety of robust, accurate, and relatively simple methods, such as the fourth-order Runge-Kutta approach. In the stochastic case, we employ a method heavily inspired by that in Ref. [84] : The method we use involves splitting up the dynamics of the stochastic PDE dy = (µy − ay 2 + D∇ 2 y)dt + σydW ,into pieces which can be solved individually. We first treat the stochastic piece,whose solution, starting from any initial configuration y(x, t), is a standard geometric Brownian motionUsing the output of the stochastic step we then use a standard Euler update rulewhere ∇ 2 y(x, t) is the appropriate discretization of the laplacian on the given lattice. For mean-field, ∇ 2 y =ȳ −y, while for a one dimensional lattice with unit spacing,Finally, we take the output of this diffusion step and at each node apply a logistic growth, according toFor a more complicated growth dynamics such as what is discussed in Section IV B, we will resort to using a routine Runge-Kutta solver, as implemented in MATLAB's ode45 function, to circumvent the lack of an available analytic solution. MATLAB implementations of all routines discussed can be found on GitHub [100].

To analyze the dynamics of the seasonal model in section III, it suffices to understand the dynamics of how the mean changes in a single growth phase. At the start of the phase, we assume a distribution of the form ρ(y) ∝ρ(y) = y −2−c D for Υ ≤ y ≤ Λ, and 0 otherwise.The evolution of each node during the growth phase follows a deterministic logistic growth. In particular, the population of a node after time T g is given bywhere y f = µ/a is the saturation for logistic growth.The average population at the end of the growth phase is then obtained asTherefore,(B3) where we define M = y f /(e µTg − 1). It should be noted that there isn't a particularly convenient closed form for the integral as presented. Moreover, since we assume 0 < c D < 1, the integral diverges as (T g ) −c D for small times, meaning that it is sometimes more convenient to split the integral asFor example, in the special case of Λ → ∞, then this manipulation would become necessary, since the second integral would take on a term of O(1) in time. This leads to an anomalous growth term, where y(t) − y(0) ∝ t c D , which emphasizes the importance of the upper cutoff.If Λ is finite, and we only care about first order in time, then we indeed can just take the expectation of both sides of the logistic growth equation and find y(T g ) − y(0) T g = µ y(0) − a y(0) 2 + O(T g ).(B4) Expanding the relevant integrals givesSince 1−c D > 0 and Λ Υ, we can trivially expand this out as