Stochastic Models for Ontogenetic Growth

1. Introduction

The ontogenetic growth pattern has been a fascinating research topic in history. Based on individual growth data sets, previous researchers scaled the time to “physiological time”to obtain ontogenetic growth curves [1,2]. Von Bertalanffy proposed a first mathematical model to describe individual growth [3,4]. Specifically, the rate of change in body mass is equal to the rate of acquisition of nutrients and energy by the growing individual minus the maintenance requirement of accumulated tissue, with appropriate parameters for converting inputs to body mass, the growth rate and maintenance rate. Von Bertalanffy assumed that tissue maintenance would be proportional to mass but that nutrient provisioning would follow a surface rule and increase with the ${\displaystyle \frac{2}{3}}$ power of mass. After this, several subsequent ontogenetic growth models have been proposed [5,6,7]. Most of these models are derived from metabolic scaling theory or allometric theory and fitted with specific species data sets. The growth rate is observed to vary considerably between species and is less variable within species [8,9,10]. West et al. proposed a general ontogenetic growth model [11]. Rather than provision nutrients with the mass power ${\displaystyle \frac{2}{3}}$, this general model assumes the scaling exponent ${\displaystyle \frac{3}{4}}$ for energy acquisition with respect to mass [12,13]. Makarieva et al. extended this model to incorporate energy conservation in ontogenetic development [14]. There are also other modifications of the general ontogenetic growth model [15]. Hou et al. developed a general model for ontogenetic growth under food restriction [16]. Shi et al. compared several ontogenetic models for crop growth [17]. Based on metabolic scaling and fitting of fish data, Sibly et al. estimated that the scaling exponent for energy acquisition is in the range $0.65-0.77$ [18]. Hatton et al. conducted a comprehensive study on scaling and body mass [19]. Recently, Escala used the heart frequency to define a time scale for ontogenetic growth that can recover a general ontogenetic growth model [20].

All these mathematical models for ontogenetic growth have useful and explanatory powers. However, these mathematical models were estimated on the basis of partial information provided by their data sets. All of them do not fully utilize the data sets in their modeling studies. We may say that the data fitting was only fitted for averages of the data sets, and so solutions or solution curves at most were an approximation of averages of the data sets. The only information from the data sets that was used for the construction of the model was the averages of the data sets. There are other data sets information that was not used. For example, the variances of the data sets were not used in all these studies. For the variance of a data set of an individual growth, the variance can be considered as a result of environmental uncertainty, environmental noise, or environmental stochastic factors. Variances may also be interpreted as variations among individuals within a species. In this study, we seek to incorporate more information provided by data sets into mathematical models for ontogenetic growth. Ito’s stochastic differential equations are a good mathematical device to handle both averages and variances when building mathematical models. There is a study that applied a stochastic differential equation for ontogenetic growth [21]. Unfortunately, it is difficult to interpret its linear noise term as data variances.

In the research of a general ontogenetic growth model [11], the authors proposed two types of growth patterns, determinate growth and indeterminate growth. However, the mathematical model itself cannot distinguish these two patterns, and there was no mathematical character that could distinguish these two growth patterns. In the present study, we construct stochastic models for ontogenetic growth in terms of stochastic differential equations. It turns out that our stochastic model is capable of distinguishing determinate growth from indeterminate growth.

The rest of the article is organized as follows. In Section 2, we recall the general model for ontogenetic growth proposed by West et al. Based on the data sets in [11], we propose two stochastic models for ontogenetic growth, one for determinate growth and one for indeterminate growth. In Section 3, we perform parameter estimation and model fitting. Firstly, we fit the average growth trajectory to a deterministic curve. Secondly, we incorporate stochasticity by identifying the best-fitted white noise term to represent growth variance. By combining these two fits, we obtain stochastic differential equation models that improve the prediction of growth trajectories, capturing the variability of ontogenetic growth. In Section 4, we provide rigorous mathematical analysis for our stochastic models. In Section 5, we make some approximations for our stochastic models and obtain a universal stochastic growth process. The article ends with a discussion section, where we propose some interesting questions and give some new perspectives on the research field.

2. Mathematical Model Setup

2.1. General Deterministic Model for Ontogenetic Growth

In this subsection, we recall the deterministic models for ontogenetic growth. We focus on a general ontogenetic growth model proposed by West et al. [11]. This model is quantitatively based on fundamental principles that at any time $t>0$, the average resting metabolic rate of an entire organism, denoted by B, can be expressed as a sum of two components:

$${\displaystyle B=\sum _c\left[N_c{B}_c+E_c{\displaystyle \frac{d{N}_c}{dt}}\right].}$$

(1)

Here, $N_c$ is the total number of cells, $B_c$ accounts for energy for the maintenance of an existing cell, and $E_c$ corresponds to the energy for the creation of a single new cell. Under the assumption that any potential variations between cells are ignored, a typical average cell mass is considered the fundamental unit $m_c$. Furthermore, it is assumed that all $E_c$, $B_c$, and $m_c$ are independent of total biomass m, and remain constant throughout growth and development processes. At any time t, we assume that $m:=m\left(t\right)=m_c{N}_c$ where $N_c=N_c\left(t\right)$ is the function of time t that represents the total number of cells of a species over time. Then equation (1) can be written as

$${\displaystyle \frac{dm}{dt}}=\left({\displaystyle \frac{m_c}{E_c}}\right)B-\left({\displaystyle \frac{B_c}{E_c}}\right)m$$

(2)

Now, if we allow the relation $B=B_0{m}^{3/4}$ to be adopted where $B_0$ is a constant for a given species, then the equation (2) can be rewritten as

$${\displaystyle \frac{dm}{dt}}=a{m}^{3/4}\left[1-{\left({\displaystyle \frac{m}{M}}\right)}^{1/4}\right],$$

(3)

where $a=B_0{m}_C/E_C$ is the combined parameter, which represents the growth rate of a species, and $m=m\left(t\right)$ is the variable that represents the total biomass of a species at time t. Supported by extensive data on mammals, birds, fish, mollusks and plants, the 3/4-power scaling, attributed to “growth spurts”, characterizes the general allometry of B from birth to reproductive maturity. The model provides the basis for deriving an allometric relationship for growth rates and the timing of life-history event. The different exponents of m in the two terms on the right-hand side of $\left(\right)$ demonstrate the mechanism that the network restricts the total number of supply units to scale differently from the total number of cells supplied. This imbalance between supply and demand ultimately limits growth, and therefore, with the growth rate a, an asymptotic maximum body size M is reached. Let $m_0$ denote the initial mass of the organism at birth. Integration of the equation results in a classical sigmoidal curve, illustrating the growth trajectory of the organism over time,

$$m\left(t\right)=M{\left\{1-\left[1-{\left({\displaystyle \frac{m_0}{M}}\right)}^{1/4}\right]e^{-at/4M^{1/4}}\right\}}^4$$

(4)

The parameters a, $m_0$, and M vary between species. In their study, they estimated these parameters for 13 species from various taxons, including mammals, birds, fish, and crustaceans, and showed that the growth trajectories of all these species align with the deterministic curves $\left(\right)$. From data analysis, the authors proposed two growth patterns: determinate growth, which attains the maximum size at a certain point in their lifespan and ceases growing, and indeterminate growth, which never reach the maximum size and exhibit a continuous growth throughout their life span. Furthermore, the solution curves seem to suggest some potential frameworks for visualizing data beyond individual species, unveiling universal growth properties across taxonomic groups.

2.2. Stochastic Models for Ontogenetic Growth

Our objective is to improve the accuracy of the model in representing growth dynamics in diverse species of animals. In practical scenarios, animal growth does not adhere solely to deterministic principles, and stochasticity plays a significant role in shaping outcomes. Factors such as genetic variation, environmental factor fluctuations, and chance events introduce randomness into growth processes. Therefore, incorporating stochastic elements into the model will enable us to create enhanced models that more faithfully reflect the complexities of nature, enhancing our understanding of animal growth dynamics, and improving our ability to predict and manage biological systems in a variety of contexts.

Now, let us start with the equation (3). To facilitate the complexity of this equation, we introduce the variable $R\left(t\right)=1-{\left(m\left(t\right)/M\right)}^{1/4}$ which stands for the relative proportion of total available metabolic energy expended on growth of a species. Then differentiating $R\left(t\right)$ and using equation (3) give us the following simple equation.

$$dR\left(t\right)=-{\displaystyle \frac{a}{4M^{1/4}}}R\left(t\right)dt.$$

(5)

Noise is well known in the process of energy transformation and the variation in metabolic energy expenditure during cell growth is a key factor influencing the diversity of growth patterns. Observe that growth parameters a and M vary among different species and also carry “noises”. Therefore, we aim to represent this variation by perturbing the ratio of these two parameters with white noise. The question is how to find a good representation for the noise intensity.

We tackle this question by examining the impact of noise through studying the significant proportion of total metabolic energy allocated to growth throughout an organism’s lifetime. To facilitate the following process, we introduce the variable $r\left(t\right)=1-R\left(t\right)={\left(m\left(t\right)/M\right)}^{1/4}$ that represents the relative proportion of the total available metabolic energy that fuels the maintenance of a species. Suppose $0<t_d$ represents the death time of a species. Whether the individual exhibits determinate or indeterminate growth, under the assumption that in initiation most of the energy is allocated to cell growth and metabolic energy for cell maintenance is ignored, the total metabolic energy $E_{tot}$ expended since birth to a time $0<t<t_d$ can be calculated by integrating B up to that time. The details are as follows:

$$\begin{array}{cc}\hfill E_{tot}& ={\int }_0^{t}B\left(s\right)ds={\int }_0^t{B}_0{m}^{3/4}\left(s\right)ds={\int }_0^t{B}_{0}r{\left(s\right)}^3{M}^{3/4}ds\hfill \\ \hfill & ={\int }_0^t{B}_{0}r{\left(s\right)}^3{M}^{3/4}·{\displaystyle \frac{4M^{1/4}}{a}}·{\displaystyle \frac{dr\left(s\right)}{1-r\left(s\right)}}={\displaystyle \frac{4B_{0}M}{a}}{\int }_0^t{\displaystyle \frac{r{\left(s\right)}^3}{1-r\left(s\right)}}dr\left(s\right)\hfill \\ \hfill & ={\displaystyle \frac{4B_{0}M}{a}}\left[-{\displaystyle \frac{1}{3}}r{\left(s\right)}^3-{\displaystyle \frac{1}{2}}r{\left(s\right)}^2-r\left(s\right)-ln\left(1-r\left(s\right)\right)\right]{|}_0^t\hfill \\ \hfill & \approx {\displaystyle \frac{4B_{0}M}{a}}\left[-{\displaystyle \frac{1}{3}}{r}^3\left(t\right)-{\displaystyle \frac{1}{2}}r{\left(t\right)}^2-r\left(t\right)-ln(1-r\left(t\right))\right].\hfill \end{array}$$

Meanwhile, the allocation of resources to support cell growth within the same time frame depends on how the rate of change of total cell counts in the organism evolves over time. This allocation can be determined by examining the following integral

$$\begin{array}{cc}\hfill E_g={\int }_0^t{E}_C\left({\displaystyle \frac{d{N}_C}{dt}}\right)dt={\int }_0^t{\displaystyle \frac{E_C}{m_C}}dm& ={\displaystyle \frac{B_0}{a}}{\displaystyle \frac{E_C}{m_C}}\left(m\left(t\right)-m_0\right)\approx {\displaystyle \frac{B_{0}M}{a}}{r}^4\left(t\right).\hfill \end{array}$$

Consequently, considering that the mass ratio remains within a range between 0 and 1 for any time $t\ge 0$, we can examine the ratio of two quantities $E_{tot}$ and $E_g$ as

$$\begin{array}{cc}\hfill {\displaystyle \frac{E_g}{E_{tot}}}& \approx {\displaystyle \frac{r^4\left(t\right)}{4\left[-\frac{1}{3}{r}^3\left(t\right)-\frac{1}{2}{r}^2\left(t\right)-r\left(t\right)-ln(1-r\left(t\right))\right]}}\hfill \\ \hfill & \le {\displaystyle \frac{1/4}{1-r\left(t\right)+min\left[-1-ln(1-r\left(t\right))-\frac{1}{3}{r}^3\left(t\right)-\frac{1}{2}{r}^2\left(t\right)\right]}}\hfill \\ \hfill & ={\displaystyle \frac{1/4}{R\left(t\right)+min\left[-1-ln\left(R\left(t\right)\right)-\frac{1}{3}{r}^3\left(t\right)-\frac{1}{2}{r}^2\left(t\right)\right]}}.\hfill \end{array}$$

(6)

For a species, the energy ratio $E_g/E_{tot}$ represents the percentage of its total energy in life that is spent on cell growth. This ratio should be different for different species. From the growth data set of 13 species in [11], we observed that there are 3 species (Salmon, Cod, and Shrimp) belonging to indeterminate growth type and the remaining 10 species belonging to determinate growth type. From the fitted results in [11], it is noticeable that the indeterminate group has a much smaller fitted growth rate a and birth biomass $m_0$. Thus, the data variation in this group is much smaller than that in determinate group. Because of this, we adopt two functions $g_1\left(R\right)$ and $g_2\left(R\right)$ to represent the intensity of noise for the determinate and indeterminate groups, respectively. The forms of these two functions are based on the ratio $E_g/E_{tot}$, since it can estimate the intensity of noise in the percentage of energy expended for cell growth during energy transformation.

For simplicity, first we approximate the right-hand side of inequality (6) with a function of the form $g_1\left(R\right)=\sigma /\left(\delta +R\right)$ for the determinate group, where $\sigma$ and $\delta$ are some constants that keep the inequality right. Therefore, we adopt the function $g_1\left(R\right)$ to represent the intensity of noise in the determinate group and introduce the following perturbation method to the growth rate ratio

$${\displaystyle \frac{-a}{4M^{1/4}}}\to {\displaystyle \frac{-a}{4M^{1/4}}}+g_1\left(R\left(t\right)\right){\displaystyle \frac{dW\left(t\right)}{dt}},$$

where $W\left(t\right)$ is a standard Wiener process. With this setup, equation $\left(\right)$ becomes a stochastic differential equation that describes the dynamics of the proportion of energy expended for cell growth in determinate group

$$dR\left(t\right)={\displaystyle \frac{-a}{4M^{1/4}}}R\left(t\right)dt+{\displaystyle \frac{\sigma R\left(t\right)}{\delta +R\left(t\right)}}dW\left(t\right).$$

(7)

Taking into account that the five parameters a, M, $m_0$, $\sigma$ and $\delta$ vary between different species in the determinate group. For the indeterminate group, since this group shows a smaller data variation compared with the determinate group, we adopt the function $g_2\left(R\right)=\sigma (1-R)/(\delta +R)$ to represent the noise intensity in the indeterminate group. Then we incorporate this function into the growth rate ratio as follows

$${\displaystyle \frac{-a}{4M^{1/4}}}\to {\displaystyle \frac{-a}{4M^{1/4}}}+g_2\left(R\left(t\right)\right){\displaystyle \frac{dW\left(t\right)}{dt}}.$$

This leads to another stochastic differential equation that describes the dynamics of the proportion of energy used for cell growth in the indeterminate group

$$dR\left(t\right)={\displaystyle \frac{-a}{4M^{1/4}}}R\left(t\right)dt+{\displaystyle \frac{\sigma R\left(t\right)\left[1-R\left(t\right)\right]}{\delta +R\left(t\right)}}dW\left(t\right).$$

(8)

In Section 4, we will analyze theoretically these two stochastic equations (7) and (8) and explain why we choose these two equations for understanding the data set of 13 species in [11] from mathematical point of view. In the next section, we will showcase how to fit the growth data set of 13 species using two proposed stochastic differential equations.

3. Parameter Fitting

3.1. Dataset

In this study, we used a diverse data set consisting of 206 samples collected from a wide range of animal species, encompassing 13 different types, including mammals, birds, fish, and crustaceans. The species represented in the data set are Guppy, Heron, Hen, Salmon, Rat, Atlantic Cod, Robin, Shrew, Guinea Pig, Domestic Rabbit, Domestic Pig, Cow and Shrimp. These species were chosen to provide a comprehensive representation of various growth patterns that include both determinate and indeterminate growth. The data set provides measurements of time (denoted as t and reported in days) and organism size (denoted as m and reported in grams). Time serves as the independent variable, representing the duration of growth observation, while the size of the organism serves as the dependent variable, representing the size of the organism at each time point.

3.2. Fitting Method

We start with the relative rate of change ${\displaystyle \frac{\dot{R}}{R}}$ in the proportion of energy expended on growth of a species. Then we allow ${\displaystyle \frac{\dot{R}}{R}}=f\left(R\right)+g\left(R\right)\xi \left(t\right)$ where $f\left(R\right)$ represents the intrinsic deterministic growth rate, $\xi \left(t\right)$ represents white noise, and $g\left(R\right)$ represents the strength of white noise. To parameterize the model, we will need to perform two data fitting steps to determine $f\left(R\right)$ and $g\left(R\right)$ from our data, separately [22]. First, we use equation (4) to fit the deterministic intrinsic growth rate of each species, where parameter values of a, $m_0$, and M vary within the taxon. Then we perform scatter plots of m versus t for different species and record the values of the three parameters. There is a significant positive correlation between life span and individual biomass. From the fitted results, $f\left(R\right)={\displaystyle \frac{-a}{4M^{1/4}}}$ with the initial value $R\left(0\right)=1-{\left(m_0/M\right)}^{1/4}$ was the best fit.

Next, we fit the strength of the white noise. We work with the mass ratio R against time t. Based on different species, we put 206 data points of form $\left(t,m\left(t\right)\right)$ into 13 groups $D_i$ for $i=1,...,13.$ For each group $D_i$, we compute the standard deviation $s_j\left(t\right)=\left|R\left(t\right)-R_{estimate}\left(t\right)\right|$. As results, we obtain 206 computed pairs of the form $\left(t,s_j\right)$. Then, we utilize the function $g_1\left(R\right)=\sigma /(\delta +R)$ to fit the computed pairs belonging to the determinate group and the function $g_2\left(R\right)=\sigma (1-R)/(\delta +R)$ to fit the computed pairs belonging to indeterminate group. The values of parameters $\sigma$ and $\delta$ are recorded for each group.

3.3. Drift Term Fitting

The data set comprises points in the form $\left(t,m\left(t\right)\right)$, and the units of time and mass are day and gram, respectively. Each group of species is examined separately. Scatter plots of biomass m versus t depicting growth patterns of selected species are analyzed. The solution curve $\left(\right)$ is used to fit the deterministic intrinsic growth rate a, the initial mass $m_0$ and the maximum biomass M. These values are listed in Table 1. Notice that there are 3 species which are of indeterminate growth types: Salmon, Cod, and Shrimp. The parameter values for these indeterminate growth species we obtained are different from those presented in the West et al. paper [11].

For 13 different animal species, we fit the data using equation (4). In Figure 1, we give four examples of fitting curves, where two are determinate growth and the other two are indeterminate growth.

Furthermore, the solution curves suggest some potential frameworks for visualizing data beyond individual species, unveiling universal growth properties across taxonomic groups. In fact, solving equation (5) gives $R\left(t\right)=R\left(0\right)exp\{-{\displaystyle \frac{-at}{4M^{1/4}}}\}$. If the dimensionless mass ratio R is plotted against a dimensionless time variable $\tau ={\displaystyle \frac{at}{4M^{1/4}}}-ln\left(R\left(0\right)\right)$, the resulting graph predicts that the body size of all species should adhere to a parameterless universal curve of $R\left(\tau \right)=e^{-\tau }$. The slopes for the growth curves of the 13 animal species closely approximate $-1$, which implies that growth decelerates as body size increases due to limitations in resource provision by biological networks. Figure 2 shows a plot of the universal deterministic growth curve.

3.4. Diffusion Term Fitting

The patterns of determinate growth species are characterized by a continuous and bounded trajectory. In contrast, the indeterminate growth group exhibits unbounded growth punctuated by intermittent interruptions during the reproduction process. The two groups of growth pattern can be fitted well with the two proposed stochastic models for ontogenetic growth.

Noise intensity

In the following data fitting process, we employ the re-scaled time $\tau =-{\displaystyle \frac{a}{4M^{1/4}}}t$ in conjunction with the metabolic distribution $R\left(\tau \right)=1-{\left({\displaystyle \frac{m\left(\tau \right)}{M}}\right)}^{1/4}$. After applying the time-rescaling technique, the stochastic differential equations $\left(\right)$ and $\left(\right)$ are transformed into the following forms.

$$dR\left(\tau \right)=-R\left(\tau \right)d\tau +{\displaystyle \frac{\overline{\sigma }R\left(\tau \right)}{\overline{\delta }+R\left(\tau \right)}}dW\left(\tau \right),$$

(9)

and

$$dR\left(\tau \right)=-R\left(\tau \right)d\tau +{\displaystyle \frac{\overline{\sigma }R\left(\tau \right)(1-R\left(\tau \right))}{\overline{\delta }+R\left(\tau \right)}}dW\left(\tau \right),$$

(10)

where $\overline{\sigma }$ and $\overline{\delta }$ are positive constant. The detailed proof of the transformation will be presented later in Section 5.

This transformation involves the conversion of data points from the study by West et al., originally represented as $\left(t,m\left(t\right)\right)$, into an alternative format of $\left(\tau ,R\left(\tau \right)\right)$. Subsequently, we proceed with a comparative analysis, aligning these data points with their corresponding estimated values which are generated from the output of the universal growth curve $R\left(\tau \right)=e^{-\tau }$. We compute the standard deviation $\left|R_{real}\left(\tau \right)-R_{estimate}\left(\tau \right)\right|$, the resulting values are then assigned labels $s_i$ ($i=1,2,...$) according to the order of the data samples provided within the dataset.

For species with determinate growth patterns, our objective is to ensure that white noise effectively captures the variability within the data set; therefore, we seek to adjust the noise intensity function $g_1\left(R\right)={\displaystyle \frac{\overline{\sigma }}{\overline{\delta }+R}}$ to align with the observed standard deviations $(\tau ,s_i)$. The optimal fit occurs when $\overline{\sigma }=0.045$ and $\overline{\delta }=0.9667$, resulting in $RMSE=0.0215$. For a group of organisms with indeterminate growth patterns, we fit the standard deviations $(\tau ,s_i)$ to the noise intensity function $g_2\left(R\right)={\displaystyle \frac{\overline{\sigma }(1-R)}{\overline{\delta }+R}}$. The optimal fit occurs when $\overline{\sigma }=0.0181$ and $\overline{\delta }=0.2843$, resulting in $RMSE=0.0426$. Figure 3 shows the fit of the dimensionless standard deviation versus the dimensionless time.

Next, a similar procedure will be performed individually for each animal group to ensure precision. One point to note is that the metabolically distributed $R\left(t\right)$ is employed but no time rescaling will be carried out; we transform the original data into the format of $\left(t,R\left(t\right)\right)$. We proceed by aligning these real data points with their corresponding estimated value which is generated from the output of the deterministic solution curve $R\left(t\right)=R\left(0\right)·exp\{-{\displaystyle \frac{at}{4M^{1/4}}}\}$. We compute the standard deviation $\left|R_{real}\left(\tau \right)-R_{estimate}\left(\tau \right)\right|$ and present the results as scatter plots. The fitted function $g_1\left(R\right)=\sigma /(\delta +R)$ is used to fit the standard deviations of the determinate group data set. Figure 4 shows four examples of standard deviation fits. These four species are all of the determinate growth type.

For the indeterminate growth, we take $g_2\left(R\right)=\sigma (1-R)/(\delta +R)$ to fit the standard deviations for the data set of the indeterminate group. Figure 5 shows two examples of the fitted results.

The optimal fits for the parameters $\sigma$ and $\delta$ for each species with determinate growth are listed in Table 2.

The optimal fits for the parameters $\sigma$ and $\delta$ for each species with indeterminate growth are listed in Table 3.

Stochastic solution paths

From the stochastic differential equation (7) for determinate growth, Ito’s formula gives us the solution as follows

$$R\left(t\right)=R\left(0\right)exp\left\{{\displaystyle -\frac{at}{4M^{1/4}}-{\int }_0^t\frac{{\sigma }^2}{2{\left(\delta +R\left(s\right)\right)}^2}ds+{\int }_0^t\frac{\sigma }{\delta +R\left(s\right)}dW\left(s\right)}\right\}.$$

(11)

Based on 206 biomass data, we use $R\left(0\right)$ as an initial condition for SDE and simulate 10 solution sample paths for each species (which makes up to 100 in total). By comparing the data points of each species with the stochastic solution paths corresponding to the set of fitted parameters $\left(a,m_0,M,\sigma ,\delta \right)$, we can see that the above equation (11) captures well the growth pattern of 10 species in the determinate growth. Equation (11) generates a band for each species, with a central curve representing the average biomass and a width indicative of the variance in mass over time.

From these solution sample paths, we obtain an empirical distribution at time t that provides the mean tumor size ${\mu }_t$ and standard deviation ${\sigma }_t$. We check if $\left|R\left(t\right)-{\mu }_t\right|\le 3{\sigma }_t$, which gives us a confidence interval of 99%. The results are shown below, which show how the real body sizes will fall in the region of our model prediction if we give a 99% confidence level.

For three species exhibiting indeterminate growth (Salmon, Cod, and Shrimp), when an organism reaches the age for reproduction, a substantial portion of its metabolic energy is allocated to reproduction, resulting in a reduction in the growth rate. We use the stochastic differential equation (8) to catch up on the growth dynamics of these species. To demonstrate the robustness of capturing the noise intensities in the data set of these species using the function $g_2\left(R\right)$, we use both stochastic differential equations (7) and (8) to fit these data to see which one gives better fitted results. Figure 7 shows the fitted solution paths of (7) and (8) for species in an indeterminate group.

4. Theoretical Analysis

In this section, we will examine the theoretical behaviors of two proposed stochastic differential equations (7) and (8) and explain why these equations are biologically relevant for fitting the data set of 13 species.

Assume that we are working on a complete probability space $\left(\Omega ,\mathcal{F},{\left\{{\mathcal{F}}_t\right\}}_{t\ge 0},\mathbb{P}\right)$ with filtration ${\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}$ satisfying the usual conditions. Since the drift term $-{\displaystyle \frac{aR}{4M^{1/4}}}$ and the two diffusion terms ${\displaystyle \frac{\sigma R}{\delta +R}}$ and ${\displaystyle \frac{\sigma R(1-R)}{\delta +R}}$ are continuously differential functions and hence locally Lipschitz continuous, there exists a unique locally almost surely continuous solution $R\left(t\right)$ for both equations (7) and (8). Using standard arguments as in [23], these two equations have unique globally almost surely continuous solutions $R\left(t\right)$ in $(0,\infty )$ whenever their initial conditions start in $(0,\infty )$. The following two theorems will summarize the dynamical behavior of the solutions of the two equations.

Theorem 1.

For any initial value $R\left(0\right)\in (0,\infty )$, there exists a unique almost surely continuous solution $R\left(t\right)$ to the equation (7) that remains in $(0,\infty )$ for all time $t\ge 0$, and $R\left(t\right)$ is a strong Markov process that satisfies the Feller property. Furthermore, if $R\left(0\right)\in (0,\infty )$, then $R\left(t\right)$ converges exponentially to 0 a.s.

Proof. By applying the Ito’s formula to the equation (7), we obtain

$$\begin{array}{cc}\hfill R\left(t\right)& =R\left(0\right)exp\left\{{\displaystyle -\frac{at}{4M^{1/4}}-{\int }_0^t\frac{{\sigma }^2}{2{\left(\delta +R\left(s\right)\right)}^2}ds+{\int }_0^t\frac{\sigma }{\delta +R\left(s\right)}dW\left(s\right)}\right\}\hfill \end{array}$$

Because the Brownian Motion $W\left(t\right)$ can take either negative values or positive values with positive probability, so without loss of generality we can assume that $\sigma >0$. Note that all the other parameters are assumed to be positive. Since $R\left(t\right)>0$ a.s. for all $t>0$, the function ${\displaystyle \frac{\sigma }{\delta +R\left(s\right)}}$ is bounded and hence, by law of large number, we get

$$\underset{t\to \infty }{lim}{\displaystyle \frac{1}{t}}{\int }_0^t{\displaystyle \frac{\sigma }{\delta +R\left(s\right)}}dW\left(s\right)=0a.s.$$

On the other hand,

$${\displaystyle \frac{1}{t}}{\int }_0^t{\displaystyle \frac{{\sigma }^2}{2{(\delta +R\left(s\right))}^2}}ds$$

is always bounded for all time $t>0$. Therefore, $R\left(t\right)$ converges exponentially to 0 a.s.

Theorem 2.

For any initial value $R\left(0\right)\in (0,1)$, there exists a unique almost surely continuous solution $R\left(t\right)$ to the equation (8) that remains in $(0,1)$ for all time $t\ge 0$, and $R\left(t\right)$ is a strong Markov process that satisfies the Feller property. Furthermore, if $R\left(0\right)\in (0,1)$, then $R\left(t\right)$ converges exponentially to 0 a.s.

Proof. The argument is totally similar to the proof of Theorem 1. Note that, when $R\left(0\right)\in (0,1)$, $R\left(t\right)\in (0,1)$ a.s. for all times t and hence the function ${\displaystyle \frac{\sigma (1-R(t\left)\right)}{\delta +R\left(t\right)}}$ is always bounded.

We will end this section with several remarks. First, the two above theorems guarantee that the solutions of both equations (7) and (8) are quickly drawn to the boundary point 0. This is an important property for both equations to fit the data set of 13 species well because if we translate the data set into the form $(t,R(t\left)\right)$ then the data will show a tendency towards zero. Biologically, this makes sense since, over time, each species uses fewer and fewer metabolic energies for its cell growth, and finally this energy decays to 0. Second, the diffusion term $g_2\left(R\right)$ in (8) differs from $g_1\left(R\right)$ in (7) by a factor $1-R$. We use this factor in $g_2\left(R\right)$ to dampen the function $g_1\left(R\right)$ since the data variation in the indeterminate growth group is much smaller than in the determinate growth group. The other reason is that it makes sure that the solution of equation (8) will live within the interval $(0,1)$ whenever it starts in $(0,1)$ and therefore it will capture the small variations in the indeterminate group.

5. Universal Stochastic Growth Process

Since growth patterns vary between different species that underlie the difference in the parameter set $(a,m_0,M,\sigma ,\delta )$, our aim is to propose a universal stochastic process that is capable of encapsulating growth patterns between diverse species. To perform this task, we employ the technique of random time change for a stochastic process.

Theorem 3

(Time change formula for Ito integral, see [26]). Suppose $c(s,\omega )$ and $\alpha (s,\omega )$ are s-continuous functions such that $\alpha (0,\omega )=0$ and $E\left[{\alpha }_t\right]<\infty .$ Let $B\left(s\right)$ be an m-dimensional Brownian motion and $v(s,\omega )$ is bounded, measurable and s-continuous. Define

$${\displaystyle \tilde{B}\left(t\right)=\underset{k\to \infty }{lim}\sum _i^k\sqrt{c({\alpha }_j,\omega )}\Delta B\left({\alpha }_j\right)={\int }_0^{{\alpha }_t}\sqrt{c(s,\omega )}dB\left(s\right).}$$

Then $\tilde{B}\left(t\right)$ is a Brownian motion and martingale with respect to ${\mathcal{F}}_{{\alpha }_t}$ and

$${\displaystyle {\int }_0^{{\alpha }_t}v(s,\omega )dB\left(s\right)={\int }_0^{t}v({\alpha }_r,\omega )\sqrt{{\alpha }^{\prime }(r,\omega )}d\tilde{B}\left(r\right)}$$

for almost surely P, where ${\alpha }^{\prime }(r,\omega )$ is derivative of $\alpha (r,\omega )$ with respect to r and so that

$${\alpha }^{\prime }(r,\omega )={\displaystyle \frac{1}{c({\alpha }_r,\omega )}},\mathrm{for}\mathrm{a}.\mathrm{a}.r\ge 0,\omega \in \Omega .$$

First, we apply the theorem to the equation (7). Note that

$$v(t,\omega )={\displaystyle \frac{\sigma R\left(t\right)}{\delta +R\left(t\right)}}$$

is measurable, bounded and continuous with respect to time t. We consider the inverse-time change

$$\alpha \left(t\right)=\tau \left(t\right):={\displaystyle \frac{at}{4M^{1/4}}}.$$

We can easily check that the time-rescale function is t-continuous, has output of 0 at initial time $t=0$, and for any fixed t, we have $E\left[\tau (t,R)\right]={\displaystyle \frac{at}{4M^{1/4}}}<\infty .$ Furthermore, by computing the derivative of $\alpha \left(t\right)$ with respect to time t, it yields the rate of change $c\left(t\right)={\displaystyle \frac{4M^{1/4}}{a}}.$ Therefore, we can define a time-changed Brownian motion

$${\displaystyle \tilde{W}\left(t\right)={\int }_0^{\tau \left(t\right)}\sqrt{c\left(s\right)}dW\left(s\right)={\int }_0^{\tau }\sqrt{{\displaystyle \frac{4M^{1/4}}{a}}}dW\left(s\right)=\sqrt{{\displaystyle \frac{4M^{1/4}}{a}}}W\left(\tau \right).}$$

Then

$${\displaystyle {\int }_0^{\tau }{\displaystyle \frac{\sigma R\left(s\right)}{\delta +R\left(s\right)}}dW\left(s\right)={\int }_0^t{\displaystyle \frac{\sigma R\left(\tau \right)}{\delta +R\left(\tau \right)}}\sqrt{{\displaystyle \frac{a}{4M^{1/4}}}}d\tilde{W}\left(s\right)={\int }_0^t{\displaystyle \frac{\sigma R\left(\tau \right)}{\delta +R\left(\tau \right)}}dW\left(\tau \right).}$$

Under the random time change, the equation (7) becomes a new SDE

$$dR\left(\tau \right)=-R\left(\tau \right)d\tau +{\displaystyle \frac{\sigma R\left(\tau \right)}{\delta +R\left(\tau \right)}}\sqrt{{\displaystyle \frac{4M^{1/4}}{a}}}dW\left(\tau \right)=-R\left(\tau \right)d\tau +{\displaystyle \frac{\overline{\sigma }R\left(\tau \right)}{\overline{\delta }+R\left(\tau \right)}}dW\left(\tau \right),$$

where $\overline{\sigma }=\sigma \sqrt{{\displaystyle \frac{4M^{1/4}}{a}}}\mathrm{and}\overline{\delta }=\delta$. By Ito’s formula, the equation (7) is equivalent to

$$R\left(\tau \right)=R\left(0\right)exp\left\{{\displaystyle -\tau -{\int }_0^{\tau }\frac{{\overline{\sigma }}^2}{2{\left(\overline{\delta }+R\left(s\right)\right)}^2}ds+{\int }_0^{\tau }\frac{\overline{\sigma }}{\overline{\delta }+R\left(s\right)}dW\left(s\right)}\right\}$$

With the same inverse-time change $\tau \left(t\right)$ and the same argument as above, the equation (8) is equivalent to

$$R\left(\tau \right)=R\left(0\right)exp\left\{{\displaystyle -\tau -{\int }_0^{\tau }\frac{{\overline{\sigma }}^2{(1-R\left(s\right))}^2}{2{\left(\overline{\delta }+R\left(s\right)\right)}^2}ds+{\int }_0^{\tau }\frac{\overline{\sigma }(1-R\left(s\right))}{\overline{\delta }+R\left(s\right)}dW\left(s\right)}\right\}.$$

It is challenging to derive a universal stochastic process from these equations after time change. Based on the theoretical result in Section 4, we know that both solutions to these two equations will converge exponentially to 0 a.s. Thus, we propose an estimation as follows. As $R\left(\tau \right)$ approaches 0 and as $\tau$ is large, we approximate the diffusion terms ${\displaystyle \frac{\overline{\sigma }R\left(\tau \right)}{\overline{\delta }+R\left(\tau \right)}}$ of (7) and ${\displaystyle \frac{\overline{\sigma }R\left(\tau \right)(1-R\left(\tau \right))}{\overline{\delta }+R\left(\tau \right)}}$ of (8) by a single term ${\displaystyle \frac{\overline{\sigma }}{\overline{\delta }}}R\left(\tau \right)$. Hence, we propose

$$dR\left(\tau \right)=-R\left(\tau \right)d\tau +KR\left(\tau \right)dW\left(\tau \right)$$

(12)

with $K=\overline{\sigma }/\overline{\delta }$. The idea is that with a large data set we can find universal noise that can be applied to all species. With the same mean and slightly larger variance, the equation offers an adequately good stochastic process to capture changes in variation of all species. By Ito’s formula, we can solve for

$$R\left(\tau \right)=R\left(0\right)exp\left\{\left(-1-{\displaystyle \frac{K^2}{2}}\right)\tau +KW\left(\tau \right)\right\}$$

(13)

Let $T=\tau +{\displaystyle \frac{K^2}{2}}\tau$, by performing another time change, equation $\left(\right)$ gives

$$R\left(T\right)=R\left(0\right)exp\left\{-T+{\displaystyle \frac{1}{\sqrt{1+\frac{K^2}{2}}}}\overline{W}\left(T\right)\right\}$$

(14)

with $\overline{W}\left(T\right)$ is a standard Brownian motion. And we can evaluate ${\displaystyle \frac{1}{\sqrt{1+\frac{K^2}{2}}}}\overline{W}\left(T\right)$ as a Brownian motion $B\left(T\right)$ with mean ${\displaystyle \frac{1}{1+\frac{K^2}{2}}}$ to get a universal stochastic process of the form

$$R\left(T\right)=exp\left\{-T+B\left(T\right)\right\}.$$

(15)

Given that two of our diffusion models lack closed-form solutions, we cannot obtain exact transition probability distributions for predicting organism sizes. However, we can make numerical predictions for sizes within the prescribed probability bounds. In theory, based on the proposed stochastic differential equations (7) and (8) together with their initial conditions, we can establish a confidence region of a specified probability for an individual’s size at any specific point in the future.

By running data fitting for 13 species for the equation, we have estimated K for each species listed in the following table $\left(\right)$. Then by taking the mean of the fitted K, we can find a value of K that presents the universal noise for the 13 species. Figure 8 shows the general stochastic band with the fitted parameter $K=-0.46$.

6. Discussion

This research was motivated by the broad idea that there should be a general stochastic model for ontogenetic growth if there was a general deterministic model for ontogenetic growth. We carefully studied most of the mathematical models for ontogenetic growth in history [6]. We decided to start with the model proposed by West et al. [11]. Based on the data sets in this study, we attempted to find the best fit for variances of the data under the condition that the fit should have two parameters and a simple mathematical form. It is surprising that this variance fit distinguishes two types of growth, determinate growth, and indeterminate growth. It turns out that we must have two variance fits, one for determinate growth and one for indeterminate growth. We do not achieve a general stochastic model for ontogenetic growth; instead, we obtain a general stochastic model for determinate ontogenetic growth and a general stochastic model for indeterminate ontogenetic growth if we could claim this, as West et al. did. However, we still can ask with what new variance form can we obtain a general stochastic model for ontogenetic growth based on the data sets in the study of West et al. This is an interesting question, and we will consider it in the future.

Based on their general model for ontogenetic growth, West et al. produced a universal growth curve. It is well known that, in the deterministic case, it is easy to obtain a universal solution curve after an appropriate time change. We plot such a universal growth curve as the dimensionless mass ratio R against dimensionless time. For two stochastic models for ontogenetic growth we obtained, we expect to get two universal stochastic growth curves, one for determinate growth and one for indeterminate growth. However, it was difficult to perform the stochastic time change without explicit solution forms. We may achieve them by numerical study. But, we made some approximation of noise terms and stochastic transformation of time, we then obtained an approximation of universal stochastic curve which is a geometric Brownian motion multiplied by an exponential time factor. The universal (deterministic and stochastic) growth curve should be a general phenomenon no matter what scaling exponent is used in modeling.

There were many studies about scaling law in ontogenetic growth. The researchers choose the scaling exponent ${\displaystyle \frac{3}{4}}$ or ${\displaystyle \frac{2}{3}}$ for energy acquisition with respect to mass, then perform data analysis to fit their models. Some research left the scaling exponent as an unknown and tried to estimate it from their data. It seems that a general study of the Bertalanffy-type ontogenetic growth equation is needed. For this, we may consider both metabolic power for maintenance and growth to be some scalings with respect to mass ${\displaystyle \frac{dm}{dt}}=a{m}^{\alpha }-b{m}^{\beta }$. In the classical work [6], it has already proposed such a general mathematical model. This general equation may provide more freedom to build a general stochastic model for ontogenetic growth. This is an interesting idea that is worth exploring further.

It may be interesting to study the ontogenetic growth of a cell, starting from the birth of a new cell until cell division. This will require details of cell biology and related data. The precise mathematical model may be established for a cell.

Tumor growth modeling has been a significant research area that provides not only new theoretic mathematical problems but also, more importantly, insight and treatment suggestions for cancer researchers. There are several methods to construct mathematical models for tumor growth, for example, mass conservation law, fluid dynamics. Energy allometry is a new way to build tumor growth models. Some studies have applied energy allometric theory to tumor growth [24,25]. Energy allometry may provide new insight into the field.