An attempt to unify some population growth models from first principles

Unificando Modelos de Dinâmica Populacional por Primeiros Princípios

Abstracts

In this work, some phenomenological growth models based only on the population information (macroscopic level) are deduced in an intuitive way. These models, for instance Verhulst, Gompertz and Bertalanffy-Richards models, are introduced in such a way that all the parameters involved have a physical interpretation. A model based on the interaction (distance dependent) between the individuals (microscopic level) is also presented. This microscopic model have some phenomenological models as particular cases. In this approach, the Verhulst model represents the situation in which all the individuals interact in the same way, regardless of the distance between them (mean field approach). Other phenomenological models are retrieved from the microscopic model according to two quantities: i) the way that the interaction decays as a function the distance between two individuals and ii) the dimension of the spatial structure formed by the individuals of the population. This microscopic model allows understanding population growth by first principles, because it predicts that some phenomenological models can be seen as a consequence of interaction at individual level. The microscopic model discussed here paves the way to finding universal patterns that are common to all types of growth, even in systems of very different nature.

Keywords
Population growth; Complex Systems; Mathematical Modelling

Neste trabalho, alguns modelos fenomenológicos de crescimento populacional são deduzidos de uma forma intuitiva. Estes modelos, por exemplo os modelos de Verhulst, Gompertz e Bertalanffy-Richards, são introduzidos de maneira que todos os parâmetros envolvidos tenham uma interpretação física. Um modelo baseado na interação (dependente da distância) entre os indivíduos (nível microscópico) também é apresentado. É mostrado que alguns modelos fenomenológicos são casos particulares deste modelo microscópico, de acordo com: i) a forma com que a interação competitiva entre indivíduos decai com a distância que os separam; e ii) a dimensão euclidiana da estrutura espacial formada pelos indivíduos da população. Nestas circunstâncias, o modelo de Verhulst, por exemplo, decorre quando todos os indivíduos interagem de forma independente da distância que os separam: uma situação tipo campo-médio. Este modelo microscópico permite compreender o crescimento populacional por primeiros princípios, uma vez que mostra que alguns modelos fenomenológicos são explicados a partir de interações no nível do indivíduo. Dessa forma, o modelo microscópico discutido aqui abre caminho para encontrar padrões universais que sejam comuns a todos os tipos de crescimento.

Palavras-chave
Crescimento Populacional; Sistemas Complexos; Modelagem Matemática

1. Introduction

The use of mathematical modeling to describe population growth behaviors has been of great success in the last decades. These models find application not only in Ecology [¹–³] (immediate application) but also in Sociology and Economy [⁴–⁶], among other disciplines. This wide spectrum of applicability of these models has motivated a quest for universal patterns that are present in different types of systems [⁶–¹³]. Here, the inductive process necessary to build population growth models is presented. Next, an attempt to unify such models by means of a generalized model is also discussed. The generalized model is built based on the interaction between the individuals of a population. This work aims to explain population growth from first principles.

An efficient mathematical model must predict the number of individuals of a population, as the time evolves, according to some ecological properties. This should also provide some explanation about the properties of the individuals which constitute the population, as the way they interact and how this affects the population as a whole.

A starting point to build a population growth model is to look at the empirical data. The data must give not only the first information about the behavior of the population, but also the test of the quality of the model. The model is classified as "good" if it agrees with the empirical evidences -, then it can be used to formulate hypotheses and to give some explanation about the system.

In this work, experimental data of a yeast population growth is used, dots in Fig. (1). The data were taken directly from the graph of the ref. [²]. Note that the data indicate that the yeast population grows quickly at the beginning (first 20 hours) and then saturates. The objective is to explain theoretically this behavior, having access only to these data. Also, one wants to answer the question: is it possible, only with this population information (macroscopic level), to predict what is happening with the individuals (microscopic level)? In other words, is it possible to infer the way a yeast cell is interacting with the others and with the environment?

Figure 1
Empirical data (dots) of the growth of a yeast population. The data were taken directly from the graph of the ref. [²]. The solid lines represent the predictions of theoretical models: Malthus, Verhulst, Gompertz, and generalized models. The used fit parameter values are presented in table (1). The Malthus model is good to describe the population in the initial stage of growth, but fails for longer time. Other models are very similar to describe the empirical data, but the generalized model presents the best fit, with larger $r^2$.

Mathematical models are addressed throughout this paper, from the simplest model to the more complete ones, always in an attempt to explain the population growth and using these yeast data as a validity test. The efficiency of these models, measured by the coefficient of determination$r^2$, are shown in table (1). This table also shows the parameters that the models use. These quantities will be better understood as the models are being presented through the text.

Table 1
Table of the fit parameter values used that best fit the model for the empirical data. The curves generated by these values are presented in Fig. (1). The generalized model, as expected, is the best model to describe the empirical data, characterized by the larger $r^2$.

Model	$q$	$N_0$	$k_0$	$K$	$r^2$
Malthus	1	20	$0.227$	$\infty$	–
Verhulst	1	15	0.302	1241	0.9865
Gompertz	$\to 0$	1	0.1824	1268	0.9715
Generalized	0.394	14	0.1820	1269	0.9874

To organize the ideas from the theoretical point of view, consider $N\left(t\right)$ as the population of yeast cells (the number or density of individuals) as a function of time $t$. The update of the population size, in a time interval $\mathrm{\Delta }t$, can be written as

(1) $$N\left(t+\mathrm{\Delta }t\right)=N\left(t\right)+\mathrm{\Delta }N\left(t\right).$$

Here, $\mathrm{\Delta }N\left(t\right)$ is a number representing the balance between the added (births) and removed (deaths) individuals in this time interval. Of course, if such numbers are known, one can predict exactly the population size in the next time interval. The problem is that, in general, those numbers are not known. A possible way to solve this problem is to extract information from theoretical models. In the following sections, two kinds of models to describe the population dynamics are presented. The first type, the phenomenological models, or macroscopic models, are those that take into account only information of the population as a whole, as Malthus, Verhulst, Gompertz and Bertalanffy-Richards models [¹⁴, ¹⁵]. The second type, the microscopic models, are those that take into account the properties of the system at the individual level [¹⁶]. This kind of model is built to achieve macroscopic phenomena as a consequence, or better, an emergent property, of the interactions between individuals. In this way, microscopic models allow describing population growth from first principles.

An alternative way to treat the problem is to use agent-based models (ABM) to simulate computationally the interactions between individuals, and then to verify the emergent properties of the system [¹⁷–¹⁹]. Studies going in this direction are the ones of Ref. [¹⁶, ²⁰–²²]. However, despite the ABM allows a better flexibility of forms at microscopic level, and allows predicting the macroscopic behavior from simulations, this methodology does not allow a complete understanding of the underlying mechanisms of the system. Although considering simple assumptions, the modeling here, first by ordinary differential equations and then deducing such equations from first principles, can be solved analytically, allowing a full understanding of the process.

The paper is organized in the following way. In section (2.), the easiest way to represent mathematically the process of growth is presented by the Malthus model. In section (3.), a small corruption of the Malthus model is done to obtain the Verhulst model. Besides being simple, it describes fairly well the yeast population growth data. In section (4.), an alternative model, known as Gompertz model, is proposed. This model also describe fairly well the yeast population. In section (5.), a model that has the Malthus, Verhulst and Gompertz models as a particular case is proposed, and because of this, it is called a generalized model. In section (6.), the generalized model is deduced from first principles, i.e., from the interaction between the individuals. This deduction allows to explain population growth (macroscopic level) from the way cells interact (microscopic level).

2. Malthus Model

In this section, the simplest way to describe the population growth is considered. By sheer simplicity, consider that the number of deaths and the number of births within the population of yeast are proportionate to: i) the number of individuals of the population and ii) the time interval. This consideration means that

(2) $$\mathrm{\text{Number of Births in}}\mathrm{\Delta }t=b\mathrm{\Delta }tN\left(t\right),$$

and that

(3) $$\mathrm{\text{Number of deaths in}}\mathrm{\Delta }t=d\mathrm{\Delta }tN\left(t\right),$$

where the parameters $b$ and $d$ are constants. This idea is reasonable because the larger the time interval is, the greater the numbers of births and deaths are. Likewise, the larger the population, greater is the number of births and deaths. This assumption of linearity is quite simplistic, but facilitates the analytic development. Moreover, for simplicity only, migration of cells is not considered.

The parameter $b$, the birth rate, can be interpreted as the average number of daughter cells generated by each cell of the population per time interval. This parameter is considered constant throughout the generations. The parameter $d$, the mortality rate, can be interpreted as the proportion of deaths in the time interval. Returning to the equation (1), one has

(4) $$N\left(t+\mathrm{\Delta }t\right)=N\left(t\right)+\left(b–d\right)\mathrm{\Delta }tN\left(t\right),$$

which is a discrete model to describe the dynamics of the population. Given an initial population and the parameters $b$ and $d$, the population size at time $t$ can be iteratively calculated. When $b>d$, i.e. the birth rate is higher than the mortality rate, then the population grows, otherwise the population decreases. The values of these parameters can be estimated from the empirical data which one wants to describe.

The discrete model (4) can be converted into a continuous equation. It is convenient because it allows analysis which are impossible in the discrete version. Taking the limit of infinitesimal time interval, one writes

(5) $$\underset{\mathrm{\Delta }t\to 0}{lim}\frac{N\left(t+\mathrm{\Delta }t\right)–N\left(t\right)}{\mathrm{\Delta }t}=k_{0}N\left(t\right),$$

where $k_0\equiv b–d$ is the population growth rate. Identifying the left part of the above equation as a derivative, one has an ordinary differential equation (ODE)

(6) $$\frac{dN}{dt}=k_{0}N.$$

The solution of this ODE is (see ^{appendix (A)} for details)

(7) $$N\left(t\right)=N_0{e}^{k_{0}t},$$

where $N_0$ is the initial population.

The simple model of Eq. (6), concerns an exponential growth of the population given by Eq. (7). This is the Malthus model, referring to the English economist, famous for his studies of populations and the phrase "The human population grows geometrically while the amount of food increases arithmetically" [²³]. The Malthus (exponential) model is also used to describe radioactive decay, among other multidisciplinary applications [²].

The problem is that the Malthus model, for been too simple, leads to an indefinite growth, with the population blowing up at a sufficiently long time. In fact, this result is inconsistent with what is observed by the data presented in Fig. (1), where the yeast population size presents a saturation value. The model fails to describe the population as a whole, however. Although the model fails for long times, it is very good to describe the population in the early stages of growth. In Fig. (1) one can see that the model fits the empirical data very well if one considers the growth in the first 20 hours. This result is quite commendable, given the simplicity of this model. This fact allows introducing an extremely important concept in modeling, which is the idea of validity region of a model. The comparison between the prediction of Malthus model and empirical data presented in Fig. (1) clearly shows that the model is good, provided it is applied to describe the initial dynamics of the population. Outside this range, the model no longer makes sense.

3. The Verhulst Model

The Malthus model was not good enough to describe the yeast population for a longer time because it does not consider the scarcity of resources. It is too simple and therefore one needs to add to it little more information to obtain better results.

As the Malthus model behaves well in the early stage of growth, one still relies on this model, but making a small change in it. Rewriting it with the introduction of a correction term in the ODE (6) yields to

(8) $$\frac{dN}{dt}=k_{0}N–\mathrm{\text{term}}.$$

This "term" should be zero (or near zero) when $N$ is sufficiently small (reaching the Malthus model in this regime) and should be maximum when the population reaches a certain level. Verhulst [²⁴, ²⁵] considered this corrective term to be proportional to $N^2∕K$, where $K$ is the carrying capacity of the population, which is the maximum size of the population that can be supported by the environment. Thus, with corrective term the Malthus ODE takes the form of the Verhulst model

(9) $$\frac{dN}{dt}=k_{0}N\left[1–\left(\frac{N}{K}\right)\right].$$

When the population size approaches the carrying capacity ($N\to K$), then $dN∕dt=$ 0, i.e. the population stops growing. This model predicts, when $t\to \infty$, that the population saturates instead of blowing up, as in the Malthus model.

The solution of this model is obtained by integrating both sides of Eq. (9) (see ^{appendix (A)} for details), resulting in

(10) $$N\left(t\right)=\frac{K}{1–\left(1–\frac{K}{N_0}\right)e^{–k_{0}t}}.$$

Note in the fit of Fig. (1) that the Verhulst model is appropriate to describe the yeast population growth. The solution (10) describes the population dynamics of yeast fairly well, both for early and late time growth. One can say that the Verhulst model (although its simplicity, using only two parameters: the growth rate and carrying capacity) captures the essence of the yeast population growth. This good description occurs even disregarding most of the details involved in the dynamics. Nevertheless supplying more details, more parameters are need to describe the situation and perhaps impedes the theoretical treatment, hardening the understanding of the phenomenon. With this simple version, considering only phenomenological parameters, it is possible to describe this population quantitatively, with an analytic solution. In fact, the model explains quantitatively that the yeast population is growing rapidly in the beginning and then saturates due to the scarcity of resources.

4. The Gompertz Model

In section (2), it was seen that the Malthus model fails because it leads to an unlimited growth of the population. In section (3.), a way (Verhulst model) was proposed to induce a saturation in the population dynamics in an attempt to address the problem of the Malthus model. An alternative idea is to induce the saturation of the population still using the Malthus model as the starting point.

The indeterminate growth of the Malthus model is due to the per capita growth rate $k_0$ remaining constant for all individuals of the population and along the generations. Consider that this growth rate is time dependent. This allows writing the model (6) as

(11) $$\frac{dN}{dt}=k\left(t\right)N.$$

To have a situation in which the population stops growing along time is interesting to consider $k\left(t\right)$ decreases with $t$. For example, an exponential decay of this growth rate, i.e. $k\left(t\right)=c^{–k_{0}t}$, where $c$ is a constant and $k_0$ plays here the role of a "half-life" growth rate, justifying the initial choice in Eq.(6). Thus one has the Gompertz model

(12) $$\frac{dN}{dt}=c{e}^{–k_{0}t}N,$$

in honor of the English economist who, in the XIX century, used it to describe human mortality [²⁶]¹. Integrating both sides of Eq. (12) (see ^{appendix (A)}), it is possible to find the solution

(13) $$N\left(t\right)=N_0{e}^{–c\left(\frac{1}{k_0}+\frac{1}{k}{e}^{–k_{0}t}\right)}.$$

Note that for $t\to \infty$ and $k_0>0$ the population converges to

(14) $$K=N_0{e}^{–\frac{c}{k_0}},$$

carrying capacity. This implies that $c∕k_0=–ln\left(K∕N_0\right)$ and the above solution can be written in terms of the parameters $K$, $N_0$ and $k_0$, not $c$. That is, the Gompertz model solution may also be written as

(15) $$N\left(t\right)=K{e}^{ln\left(\frac{N_0}{K}\right)e^{–k_{0}t}}.$$

A more convenient way to write the Gompertz model (12) is by the ODE

(16) $$\frac{dN}{dt}=–k_{0}Nln\left(\frac{N}{K}\right),$$

which is more usually found in the literature [²⁷–³⁰]. Note that integrating the Eq. (16) leads to Eq. (15), which shows that the Eq. (12) and Eq. (16) describe exactly the same model.

Applying the Gompertz model to describe the yeast population growth, one gets a result as good as the Verhulst model, as shown in Fig. (1). Therefore this model is also interesting to model the population dynamics of these microorganisms.

5. The Generalized Model and the Bertalanffy-Richards Model

Both the Verhulst model (which considers a saturation term) and the Gompertz model (which considers a growth rate that decreases exponentially as a function of time) lead to satisfactory results to describe the yeast population growth. But if these two essentially different models lead to quite similar results, which of them is true? Or are both false? To answer these questions, it is interesting to work out a more complete model, which has at least the Verhulst and Gompertz models as special cases.

Consider that the logarithm function of the Gompertz model (16) can be interpreted as a special case of a more general function. One option is to look at the generalized logarithm function of statistical mechanics (which considers $q+1$). This function is defined by (see ^{appendix (B)} for details)

(17) $${ln}_q\left(x\right)\equiv \frac{x^q–1}{q},$$

where $q$ is the generalization parameter. This function is called this way because it recovers the usual logarithm function when one takes the limit $q\to 0$ (see ^{appendix (B)}) ². Replacing the usual logarithm function of the Gompertz model Eq. (16) by the generalized logarithm function of Eq. (17), one writes

(18) $$\frac{dN}{dt}=–k_{0}N{ln}_q\left(\frac{N}{K}\right),$$

where $q$ plays the role of a generalization parameter. This model will be called generalized model because it retrieves not only the Gompertz model in the limit $q\to 0$, but also the Verhulst model when $q=1$. Moreover, if $q=1$ (Verhulst) and $K\to \infty$, the generalized model also recovers the Malthus model. For more details, see Ref. [¹⁴, ¹⁶, ³¹–³⁴].

The solution of the generalized model can be obtained integrating both sides of (18) (see ^{appendix (A)}), which leads to

(19) $$N\left(t\right)=\frac{K}{{\left[1–\left[1–{\left(\frac{K}{N_0}\right)}^q\right]e^{–k_{0}t}\right]}^{\frac{1}{q}}}.$$

Note that for $q=1$ and $q\to 0$, this solution leads to Eqs. (10) and (15), respectively.

This model was originally proposed in Ref. [³⁵] and represents an important step in building unified population growth models. It has been applied to the yeast experimental data, as can be seen in Fig. (1). The best fit has been obtained with $q=0.394$, which means an intermediate situation between the Gompertz ($q\to$ 0) and the Verhulst ($q=1$) models, is better to adjust these data.

The parameter $q$ was introduced here as a theoretical argument of generalization, without presenting any physical interpretation. With a microscopic approach that takes into account the interactions between individuals, see Section (6), it is possible to give a physical interpretation to the generalization parameter.

Besides Eq. (18) to be based on the macroscopic properties of the system, a phenomenological approach, it can give some hints about the average behavior of individuals. For instance, one can compute the per capita growth rate$G$, from Eq. (18), by

(20) $$G\left(N\right)=\frac{1}{N}\frac{dN}{dt}=–k_0{ln}_q\left(\frac{N}{K}\right).$$

If Eq. (19) is putted in the equation above, the time evolution of the per capita growth rate $G\left(t\right)$ is obtained. The plot of Fig. (2) presents the time evolution of this quantity. With the exception of Malthus model, all the models covered by the generalized model present a per capita growth rate which decreases as a function of time. The Verhulst, Gompertz and other particular cases differ from one another just with respect to the way that the individuals of the population decrease the capability to reproduce as time evolves.

Figure 2
Plot of the per capita growth rate as a function of time. Except the Malthus model, in which the per capita growth rate is constant throughout the generations, the other particular cases of the generalized models (as Verhulst and Gompertz models) differ only in the way that this average individual growth rate decreases as a function of time.

The per capita growth rate $G$ is an individual property which has been inferred from macroscopic information. This procedure to infer system characteristics is known as the top down approach. This quantity is not obtained by "first principles" (bottom up approach). In the next section, one introduces a microscopic model for the interaction between the cells of the population.

The generalized model, when written in a more explicit form, i.e. using Eq. (17), becomes

(21) $$\frac{dN}{dt}=–\left(\frac{k_0}{q{K}^q}\right)N^{q+1}+\frac{k_0}{q}N,$$

or, in a more compact form

(22) $$\frac{dN}{dt}=a{N}^{\beta }–bN,$$

where the parameters $a\equiv –k_0∕\left(q{K}^q\right)$, $b\equiv –k_0∕q$, and $\beta \equiv q+1$ is introduced. Eq. (22) is known in the literature as the Bertalanffy-Richards model [¹⁵,³⁵–³⁷]. This model is often used to describe ontogenetic growth [¹²] and has recently been applied to describe human population growth [⁶]. Interestingly in this identification is that it provides interpretation of the parameters of the Bertalanffy-Richards model, i.e.,

• growth rate: $k_0=–b\left(\beta –1\right)$;

• carrying capacity: $K={\left(\frac{b}{a}\right)}^{1∕\left(\beta –1\right)}$;

• generalization constant: $q=\beta –1$.

As Eq. (18) retrieves some phenomenological population growth models as particular cases, the Bertalanffy-Richards model also does that. For instance: $\beta =2$ ($q=1$) leads to Verhulst Model if the carrying capacity is finite, or to the Malthus model otherwise; and $\beta \to 1$ ($q\to 0$) leads to the Gompertz model.

However, as there is no physical interpretation of the parameter $q$ or $\beta$ yet. This interpretation are achieved in the next section with the consideration of a microscopic model based on the interactions of the individuals.

6. Microscopic Model

So far, the presented models have been built considering only information in the macroscopic level, i.e. population information. However, to address population dynamics from a deeper point of view, it is imperative to model the interactions between individuals. This microscopic (or bottom up) modeling leads to a description of population dynamics from first principles.

The first step towards this microscope modeling was presented by Mombach et al in Ref. [³⁸] and reworked in [¹⁰, ¹⁶, ³⁹]. The result of this model is the fact that some observed macroscopic properties are consequences of the inhibitory interactions between the microcomponents of the system. Yet, a physical interpretation of the generalization parameter $q$ is given in Ref. [¹⁴]. In this section, this model is discussed in detail.

To begin with, consider that the replication rate of a cell is regulated by two factors: one is intrinsic to the cell and the other is related to the interactions of the cell with the other cells of the population. This means that the following scheme for a certain cell must be valid:

[Replication rate] = [self-replicate stimulus] +

+ [stimulus from interaction with other cells].

Following this scheme and naming $R_i$ the replication rate of the $i$-th cell, the number of daughter cells generated by $i$ in a time interval $\mathrm{\Delta }t$, is

(23) $$\mathrm{\Delta }t{R}_i=G_i+J{I}_i.$$

Here, $G_i$ is the intrinsic ability of the cell to self-replicate; $I_i$ is the interaction field felt by this cell and caused by the other cells of the population and $J$ is the intensity of this interaction field. If $J>0$, there is cooperation among cells and if $J<0$, there is competition among them. In this paper, only the competition is explored. For readers interested in cooperation between cells, see Ref. [¹⁶, ³³].

The total number of daughter cells added to the population in a certain generation $t$ is $\mathrm{\Delta }N\left(t\right)=\mathrm{\Delta }t{\sum }_{i=1}^N{R}_i$, which leads to

(24) $$N\left(t+\mathrm{\Delta }t\right)=N\left(t\right)+\mathrm{\Delta }t\sum _{i=1}^N{R}_i.$$

Considering infinitesimal time intervals $\mathrm{\Delta }t\to 0$, the recurrence equation above, together with Eq. (23), becomes

(25) $$\frac{dN}{dt}=N⟨G⟩+J\sum _{i=1}^N{I}_i,$$

where $⟨G⟩\equiv \left(1∕N\right){\sum }_{i=1}^N{G}_i$ is the average value of the intrinsic replication capacity of a single cell.

The interaction field $I_i$ is the result of the interaction of a cell with all the others. This allows writing the summation $I_i={\sum }_{j\ne i}{I}_{ij}$, where $I_{ij}$ is the interaction field between cells $i$ and $j$. Consider, by hypothesis, that the interaction $I_{ij}$ decays as function of the distance $r_{ij}$ between them according to

(26) $$I_{ij}\left(r_{ij}\right)=\frac{1}{r_{ij}^{\gamma }},$$

where $\gamma$ is the decay exponent of the interaction field. For $\gamma =0$, the interaction field does not depend on the distance, i.e. the region of interaction between two cells is infinite. This case is called mean field situation. For $\gamma \to \infty$, the region of interaction vanishes, which means that cells do not interact at all. To preserve the internal structure of cells, consider $r_{ij}\ge 2r_0$, where $2r_0$ is the diameter of the cell. The drawing of Fig. (3) illustrates the interaction field decay as a function of the distance, for different values of the exponent decay. For a more generic approach about the interaction function of Eq. (26), see Ref. [³³].

Figure 3
Graph representing the interaction field between two individuals ($i$ and $j$) depending on the distance $r_{ij}$ between them, according to some decay exponent values. The quantity $r_0=1∕2$ is fixed. The larger $\gamma$ is, the faster the decay of the field with the distance is. For $\gamma =0$, the interaction region is infinite, i.e. the intensity of the interaction does not depend on the cell separation (mean field situation). For $\gamma \to \infty$, the interaction region vanishes, which means that cells do not interact.

With respect to the cell spatial distribution, consider that they form a structure with hypervolume $V_D$ in $D$ dimensions. For $D=3$, the hypervolume is the usual Euclidean volume; for $D=2$ the population is distributed forming a surface and the hypervolume is the area of this surface and for $D=1$ the population is distributed forming a straight segment and, in this case, the hypervolume is the length of this segment. Consider also $\rho \left(r\right)$ as the density of cells in the hypervolume element $d^{D}r$, positioned in $r$. This position vector has its origin in cell $i$, as depicted in Fig. (4). The number of cells in this hypervolume element is $dN=\rho \left(r\right)d^{D}r$ and if the population is continuously distributed in the space, the field $I_i={\sum }_{j\ne i}{I}_{ij}$ can be written as the integral $I_i={\int }_{V_D}{r}^{–\gamma }dN$, which becomes

(27) $$I_i={\int }_{V_D}\frac{\rho \left(r\right)}{r^{\gamma }}{d}^{D}r.$$

Figure 4
Hypervolume $V_D$ in $D$ dimensions, formed by the spatial cell (dots) distribution. Here a small group of cells is in a hypervolume element $d^{D}r$, localized by the position vector $r$, originated in the $i$-th cell.

Assuming, for simplicity only, that the population is homogeneously distributed in space, i.e. $\rho \left(r\right)={\rho }_0=\mathrm{\text{cte}}$ and using hyperspherical coordinates (see ^{appendix (C)} for details), the total field felt by this cell is

(28) $$I_i=I\left(N\right)=\frac{{\tilde{\rho }}_D}{D}{ln}_{1–\gamma ∕D}\left[\frac{N}{\left({\tilde{\rho }}_D∕D\right)}\right].$$

where ${\tilde{\rho }}_D$ is a parameter related to the density of cells and which depends only on the dimension $D$. Eq. (28) concerns the interaction field of the specific cell $i$, however the result is identical to all the other cells, regardless of their position. This is because the right side of the Eq. (28) does not depend on the "index" $i$. The interaction field in a single cell depends only on the size of the population $N$, which justifies to write $I_i=I\left(N\right)$ in Eq. (28).

From Eq. (25), one obtains

(29) $$\begin{array}{rcll}\frac{dN}{dt}& =& \left[\frac{J}{\left(1–\frac{\gamma }{D}\right)}{\left(\frac{{\tilde{\rho }}_D}{D}\right)}^{\frac{\gamma }{D}}\right]N^{2–\frac{\gamma }{D}}& \\ & +& \left[⟨G⟩–\frac{J{\tilde{\rho }}_D}{D–\gamma }\right]N\end{array}$$

which is nothing more than the Bertalanffy-Richards growth model given by Eq. (22) (and consequently the generalized model (18)) with

(30) $$a=\frac{J}{\left(1–\frac{\gamma }{D}\right)}{\left(\frac{{\tilde{\rho }}_D}{D}\right)}^{\frac{\gamma }{D}},$$

(31) $$b=\frac{J{\tilde{\rho }}_D}{D–\gamma }–⟨G⟩,$$

and

(32) $$\beta =2–\frac{\gamma }{D}.$$

Thus, one has an interpretation for the growth rate $k_0$, the carrying capacity $K$, and the generalization parameter $q$ (or $\beta$) from this microscopic model. Actually, from the above relations, one gets, from a microscopic point of view, the following interpretation for these quantities (given that there is competition between cells, that is $J<0$):

• generalization constant:

  (33) $$q=1–\frac{\gamma }{D};$$

• growth rate:

  (34) $$k_0=⟨G⟩q–\frac{J{\tilde{\rho }}_D}{D};$$

• carrying capacity:

  (35) $$K={\left(\frac{\frac{{\tilde{\rho }}_{D}J}{D}–\left(1–\frac{\gamma }{D}\right)⟨G⟩}{{\left(\frac{{\tilde{\rho }}_D}{D}\right)}^{\frac{\gamma }{D}}J}\right)}^{\frac{1}{1–\frac{\gamma }{D}}}.$$

The microscopic model allows stating that the generalization parameter $q$ is a relation between the repulsive potential range among the cells and the dimension of the structure formed by the population. This result proves to be important in the sense that it has only macroscopic quantities, although it has been obtained from microscopic premises. In other words, the global behavior emerges from local interactions between the cells which compose the population. Most importantly, this microscopic model explains, from first principles, all the phenomenological models discussed previously and gives a physical interpretation for all quantities considered.

As presented, the generalized model (18) has as particular cases the Gompertz and the Verhulst models. In the microscopic model, when the dynamic interaction takes place within a range $\gamma =0$, which means that the intensity of interaction among the cells does not depend on the distance, one has Verhulst dynamics. According to Eq. (33), one has, in this case, $q=1$. Thus, the Verhulst model can be interpreted as a mean field approach to the population dynamics. In this limit case, from Eqs. (34) and (35) one has the growth rate $k_0=⟨G⟩–J{\tilde{\rho }}_D∕D$ and carrying capacity $K=–{\tilde{\rho }}_D⟨G⟩∕\left(JD\right)$, respectively.

The Gompertz model emerges for $\gamma \to D$ ($q\to 0$). This phenomenological model takes place when the interaction among individuals decays with a range equal to the dimension of the system. Taking the limit $\gamma \to D$ in Eqs. (34) and (35), the population must grow with a rate $k_0=–J{\tilde{\rho }}_D∕D$, until it reaches the carrying capacity $K=\left({\tilde{\rho }}_D∕D\right)e^{–\frac{D⟨G⟩}{J{\tilde{\rho }}_D}}$.

The fact that the microscopic model leads to the generalized model shows that it is possible to find universal patterns in population growth - universal here in the sense that properties that are present in all kinds of growth should exist, even in essentially different systems. According to the model presented, a subtle difference in the way that the individuals interact leads to a huge difference in the ecological properties. This fact can be seen as a motivation to interpret all diversity of growth in nature as a result of a single and universal rule which describes the way of interaction of the components of the system.

7. Conclusion

Here, some important phenomenological population growth models have been discussed, with physical interpretation of all the parameters considered. The models were compared to empirical data (yeast population growth) to test their validity. It was also shown that these models can be retrieved as particular cases of a microscopic model, which takes into account the interaction (depending on the distance) between the individuals of the population. From this approach, it was verified that the Verhulst model represents a situation in which each individual of the population interacts with all the others in the same way, regardless the distance that separates them. That means the Verhulst model is a mean field interaction model. Other models, as Gompertz and the generalized one, are retrieved according to the relation between the interaction decay exponent and the dimension of the structure formed by the population.

The phenomenological models presented here can be seen as a consequence of the common individual level interaction. In this sense, this microscopic model takes an important step towards a more profound understanding of the population growth and also connects many types of growth in a single approach. Thus, it paves the way to find universal patterns, common in all types of growth, even in systems of very different nature.

Appendix A. Complete Solution of the Generalized Model

The generalized model is

(A-1) $$\frac{dN}{dt}=–k_{0}N{ln}_q\left(\frac{N}{K}\right).$$

This model is called so because it retrives the Malthus ($q=1$, $K\to \infty$), Verhulst ($q=1$) and Gompertz ($q\to 0$) models as particular cases. The solution of this model is obtained integrating both sides of above equation, that is

(A-2) $${\int }_{N_0}^N\frac{d{N}^{\prime }}{–\frac{k_0}{q{K}^q}{N}^{q+1}+\frac{k_0}{K}{N}^{\prime }}={\int }_0^{t}d{t}^{\prime },$$

where the explicit form of the generalized logarithm (Eq. (17)) is used. Solving the integral above leads to

(A-3) $$–\left(1+q\right)ln\left(N_0\right)+ln\left[e^{\left(1+q\right)ln\left(N_0\right)}–K^q{N}_0\right]=k_{0}t,$$

and

(A-4) $$\begin{array}{rcll}& & ln\left[\frac{N^{1+q}}{e^{\left(1+q\right)lnN}–K^{q}N}\right]& \\ & & –ln\left[\frac{N_0^{1+q}}{e^{\left(1+q\right)ln{N}_0}–K^q{N}_0}\right]=k_{0}t.\end{array}$$

Identifying $e^{\left(1+q\right)ln\left(N\right)}=e^{ln{N}^{\left(1+q\right)}}=N^{1+q}$, one obtains

(A-5) $$\frac{N^{1+q}}{N^{1+q}–K^{q}N}=\frac{N_0^{1+q}}{N_0^{1+q}–K^q{N}_0}{e}^{k_{0}t}.$$

Then, after some algebraic manipulation, the solution of the generalized model is

(A-6) $$N\left(t\right)=\frac{K}{{\left[1–\left[1–{\left(\frac{K}{N_0}\right)}^q\right]e^{–k_{0}t}\right]}^{\frac{1}{q}}}.$$

For $q=1$ and $K\to \infty$ (Malthus), the solution above leads to Eq. (7). For $q=1$ (Verhulst), this solution leads to Eq. (10). Finally, for $q\to 0$ (Gompertz), and using the property that $ln\left(x\right)=\underset{q\to 0}{lim}\left(x^q–1\right)∕q$, the solution leads to Eq. (15).

B. Generalized Logarithm and Exponential Functions

The logarithm function$ln\left(x\right)$ can be seen, among other interpretations, as the area below the hyperbole function $f\left(t\right)=1∕t$, given by the integration

(B-1) $$ln\left(x\right)={\int }_1^x\frac{1}{t}dt.$$

We can use this idea to build a generalized logarithm function, so that the usual logarithm function is only a particular case [⁴⁰]. Defining such generalized logarithm function as the area below the asymmetric hyperbole $f_q\left(t\right)=1∕t^{1–q}$ in the interval $t\in \left[1,x\right]$, shows that

$${\text{In}}_q(x)={\displaystyle {\int }_1^x\frac{dt}{t^1-q}=\{\begin{array}{c}\frac{x^q-1}{q}\text{for}q\ne 0\\ \text{In}(x)\text{for}q\to 0.\end{array}}$$

Thus, this function is a generalization, by introducing the parameter $q$, of the natural logarithm function. After all, it retrieves the logarithm function in the limit $q\to 0$, that is

(B-2) $$\underset{q\to 0}{lim}{ln}_q\left(x\right)=ln\left(x\right).$$

The inverse function of the generalized logarithm function is the generalized exponential function, defined by

(B-3) $$e_q\left(x\right)=\left\{\begin{array}{cc}\underset{q^{{}^{\prime }}\to q}{lim}{\left(1+q^{{}^{\prime }}x\right)}^{1∕q\prime }\hfill & \mathrm{\text{if}}qx>–1\mathrm{\text{;}}\hfill \\ \hfill & \hfill \\ 0\hfill & \mathrm{\text{otherwise.}}\hfill \end{array}\right.$$

For $q=0$, one recovers the exponential function: $e_0\left(x\right)=e^x$.

These generalized functions have shown to be important because they allow an easy handling of algebraic expressions, besides recovering particular cases [⁸, ¹⁶].

C. Calculation of $I_i$

Here, the calculation of the field felt by a single individual is performed. Consider that (hyperspherical coordinates)

(C-1) $$d^{D}r=r^{D–1}drd{\mathrm{\Omega }}_D,$$

which yields Eq. (27) to

(C-2) $$I_i={\tilde{\rho }}_D{\int }_{2r_0}^{R_{max}}{r}^{–\gamma +D–1}dr,$$

where ${\tilde{\rho }}_D\equiv {\rho }_0\int d{\mathrm{\Omega }}_D$ is a parameter related to the cell density which depends on the dimension $D$, and $\int d{\mathrm{\Omega }}_D$ is the solid angle. For $D=1$, $\int d{\mathrm{\Omega }}_1=2$; for $D=2$, $\int d{\mathrm{\Omega }}_2=2\pi$; and for $D=3$, $\int d{\mathrm{\Omega }}_3=4\pi$. The quantity $R_{max}$ is the maximum distance between two cells. To solve the integral (C-2) leads to

(C-3) $$I_i=\frac{{\tilde{\rho }}_D}{D–\gamma }\left[R_{max}^{D–\gamma }–{\left(2r_0\right)}^{D–\gamma }\right].$$

Considering for sheer convenience, and without loss of generality, that $r_0=1∕2$, the equation above is simplified to

(C-4) $$I_i=\frac{{\tilde{\rho }}_D}{D–\gamma }\left[R_{max}^{D–\gamma }–1\right].$$

The distance $R_{max}$ can be written in terms of the number of cells $N$. Note that $N={\int }_{V_D}\rho \left(r\right)d^{D}r$, which, given the above considerations, leads to

(C-5) $$N={\tilde{\rho }}_D{\int }_0^{R_{max}}{r}^{D–1}dr$$

and consequently

(C-6) $$R_{max}={\left(\frac{DN}{{\tilde{\rho }}_D}\right)}^{\frac{1}{D}}.$$

Introducing such results into (C-4), one calculates

(C-7) $$I_i=I\left(N\right)=\frac{{\tilde{\rho }}_D}{D\left(1–\frac{\gamma }{D}\right)}\left[{\left(\frac{N}{\left({\tilde{\rho }}_D∕D\right)}\right)}^{1–\frac{\gamma }{D}}–1\right]$$

or, in a more compact form, using the generalized logarithm of Eq. (17),

(C-8) $$I_i=I\left(N\right)=\frac{{\tilde{\rho }}_D}{D}{ln}_{1–\frac{\gamma }{D}}\left[\frac{N}{\left({\tilde{\rho }}_D∕D\right)}\right].$$

This is the interaction field of the specific cell $i$, however the result proves to be identical to all the other cells, regardless of their position. It is because the right side of the above equation does not depend on the "index" $i$. In fact, the interaction field in a single individual depends only on the size $N$ of the population.

Acknowledgements

We would like to acknowledge the useful and stimulating discussions with Alexandre Souto Martinez and Angelica Sousa da Mata.

References

• ^[1]
  J.D. Murray, Mathematical Biology: I. An Introduction. Interdisciplinary Applied Mathematics - Third Edition (Springer, New York, 2002).
• ^[2]
  L. Edelstein-Keshet, Mathematical Models in Biology (Society for Industrial and Applied Mathematics, Philadelphia, 2005), 1ª ed.
• ^[3]
  F. Ribeiro and N. Caticha, J. Theor. Biol. 258, 465 (2009).
• ^[4]
  M. Ausloos and B. Jardini, Advances in Complex Systems 16, 1250088 (2013).
• ^[5]
  D. Strzalka, Acta Physica Polonica B 40, 41 (2009).
• ^[6]
  L. Bettencourt, J. Lobo, D. Helbing, C. Kühnert and G.B. West, Proceedings of the National Academy of Sciences of the United States of America 104, 7301 (2007).
• ^[7]
  S. Solomon, https://arxiv.org/abs/cond-mat/9901250 (1999).
  » https://arxiv.org/abs/cond-mat/9901250
• ^[8]
  B.C.T. Cabella, A.S. Martinez and F. Ribeiro, Physical Review E 83, 061902 (2011).
• ^[9]
  M. Chester, Open Journal of Ecology 1, 77 (2011).
• ^[10]
  F.L. Ribeiro and K.N. Ribeiro, Physica A: Statistical Mechanics and its Applications 434, 201 (2015).
• ^[11]
  C. Guiot, P.G. Degiorgis, P.P. Delsanto, P. Gabriele and T.S. Deisboeck, Journal of Theoretical Biology 225, 147 (2003).
• ^[12]
  G.B. West, J.H. Brown and B.J. Enquist, Nature 413, 628 (2001).
• ^[13]
  D. Strzalka and F. Grabowski, Physica A: Statistical Mechanics and its Applications 387, 2511 (2008).
• ^[14]
  A.S. Martinez, R.S. González and C.A.S. Terçariol, Physica A: Statistical Mechanics and Its Applications, 387, 5679 (2008).
• ^[15]
  L. Yuancai, C. Marques and F. Macedo, Forest Ecology and Management 96, 283 (1997).
• ^[16]
  F.L. Ribeiro, Bulletin of Mathematical Biology 77, 409 (2015).
• ^[17]
  N. Boccara, Modeling Complex Systems (Springer-Verlag, New York, 2004).
• ^[18]
  S.F. Railsback and V. Grimm, Agent-Based and Individual-Based Modeling: A Practical Introduction (Princeton University Press, Princeton, 2011).
• ^[19]
  M. Mitchell, Complexity: A Guided Tour (Oxford University Press, Oxford, 2009).
• ^[20]
  M.A. Tsyganov, I.B. Kresteva, G.V. Aslanidi, K.B. Aslanidi, A.A. Deev and G.R. Ivanitsky, Journal of Biological Physics 25, 165 (1999).
• ^[21]
  D. Drasdo and S. Hohme, Mathematical And Computer Modelling, 7177, 1163 (2003).
• ^[22]
  A.H.O. Wada, T. Tome and M.J. Oliveira, J. Stat. Mech. 2015, P04014 (2015).
• ^[23]
  T. Malthus, An Essay on the Principle of Population as it Affects the Future Improvement of Society (J. Johnson, London, 1798), available in http://www.esp.org/books/malthus/population/malthus.pdf
  » http://www.esp.org/books/malthus/population/malthus.pdf
• ^[24]
  P. Verhulst, Nouveaux Memoires de l’Academie Royale des Sciences et Belles Lettres de Bruxelles 18, 1 (1845).
• ^[25]
  P. Verhulst, Nouveaux memoires de l’Academie Royale des Sciences et Belles Lettres de Bruxelles, 20, 1 (1847).
• ^[26]
  R. Gompertz, Philos Trans R Soc Lond 115, 513 (1825).
• ^[27]
  M. Sinai, Cancer Research 48, 7067 (1988).
• ^[28]
  C. Song and O. Kuznetsova, in: SAS Conference Proceedings: PharmaSUG, Boston, p. 20-23, 2001.
• ^[29]
  J.L. Haybittle, International journal of epidemiology 27, 885 (1998).
• ^[30]
  M. Molski and J. Konarski, Bio Systems 92, 245 (2008).
• ^[31]
  B.C.T. Cabella, F. Ribeiro and A.S. Martinez, Physica A: Statistical Mechanics and its Applications 391, 1281 (2012).
• ^[32]
  F. Ribeiro, B. Caetano and A.S. Martinez, Theory in Biosciences 133, 135 (2014).
• ^[33]
  R.V. Dos Santos, F.L. Ribeiro and A.S. Martinez, Journal of theoretical Biology 385, 143 (2015).
• ^[34]
  L.S. Dos Santos, B.C.T. Cabella and A.S. Martinez, Theory Biosci. 133, 117 (2014).
• ^[35]
  F.J. Richards, J Exp Bot 10, 290 (1959).
• ^[36]
  L. von Bertalanffy, Q. Rev. Biol. 32, 217 (1949).
• ^[37]
  F. Kozusko and M. Bourdeau, Cell Prolif. 40, 824 (2007).
• ^[38]
  J.C.M. Mombach, N. Lemke, B.E.J. Bodmann and M.A.P. Idiart, Europhysics Letters 59, 923 (2002).
• ^[39]
  A. D’Onofrio, Chaos, Solitons & Fractals 41, 875 (2009).
• ^[40]
  T.J. Arruda, R.S. González, C.A.S. Terçariol and A.S. Martinez, Phys. Lett. A 372, 2578 (2008).

• 1
  The life insurance started to become possible because of the Gompertz ideas in 1825. He realized that the probability of an adult to die in the next year increases exponentially with age. This model is still used by insurance agencies.
• 2
  Note that ${ln}_q\left(x\right)$ in Eq. (17)does not means "logarithm of $x$ in the base $q$ "; in this case, the notation ${log}_q\left(x\right)$ is employed.

Publication Dates

• Publication in this collection
  2017

History

• Received
  28 May 2016
• Reviewed
  07 Sept 2016
• Accepted
  13 Sept 2016