Authors: Ryo Oizumi, Hisashi Inaba
Categories: Research Article, Biology and Life Sciences, Engineering and Technology, Computer and Information Sciences, Physical Sciences, Research and Analysis Methods, Ecology and Environmental Sciences
Source: PLoS ONE
Various definitions of fitness are essentially based on the number of descendants of an allele or a phenotype after a sufficiently long time. However, these different definitions do not explicate the continuous evolution of life histories. Herein, we focus on the eigenfunction of an age-structured population model as fitness. The function generates an equation, called the Hamilton–Jacobi–Bellman equation, that achieves adaptive control of life history in terms of both the presence and absence of the density effect. Further, we introduce a perturbation method that applies the solution of this equation to the long-term logarithmic growth rate of a stochastic structured population model. We adopt this method to realize the adaptive control of heterogeneity for an optimal foraging problem in a variable environment as the analyzable example. The result indicates that the eigenfunction is involved in adaptive strategies under all the environments listed herein. Thus, we aim to systematize adaptive life histories in the presence of density effects and variable environments using the proposed objective function as a universal fitness candidate.
Since the publication of The Origin of Species by Charles Darwin, many biologists have believed that evolution is promoted by mutation and adaptation. Mutation is a well-known phenomenon that has been extensively studied at the molecular level. Similarly, adaptation is a widely accepted idea, and its degree is estimated by an index called “fitness,” which has been defined in several ways essentially based on the number of descendants of an allele or a phenotype after a sufficiently long time. If an allele or an individual with a mutation has greater fitness than other alleles/individuals without a mutation, the mutation will eventually dominate the population. However, fitness is not observed easily in nature; therefore, we must rely on indirect indices to analyze evolution.
Because it cannot be easily observed in nature, fitness does not have a unique and quantitative definition. An adaptive gene must meet several requirements to thrive in a population. The indicator must be a measure by which adaptive genes dominate the population, regardless of population dynamics, including saturated growth, exponential growth, or stochastic growth.
Biologists use population growth rates, basic reproductive numbers, and abstract payoff functions instead of substantial fitness and often call them “fitness.” These indices can represent fitness in restricted environments; e.g., (1) absence of intraspecific and interspecific competition, (2) population dynamics limited to one generation, and (3) negligible population dynamics. However, these conditions are unusual in the natural world. Therefore, the environments surrounding organisms are believed to comprise combinations of these conditions.
For example, for the combination of (1) and (2), we can determine the fitness associated with the life schedules of individuals and the population dynamics. A study addressing this problem linked age-structured models to control theory [1]. The researchers used the characteristic function of the Euler–Lotka equation as the fitness metric. Although this model was constructed to maximize the characteristic function with the adaptive life schedule, it maximized the population growth rate. The maximization of the characteristic function is equivalent to the maximization of population growth. Further, the model provided a framework for the analysis of the adaptive control of life history to natural selection.
The systematization of mathematical models related to the evolution of life histories has been promoted by linking the behavior of individuals to their population dynamics. One of the most challenging aspects of finding a general definition of fitness is that general population dynamics contain intra- and inter-specific competition, which complicates the dynamics and makes it challenging to identify what the species optimizes. The r/K selection theory argues that the adaptive life history in a stationary population maximizes the carrying capacity [2]. Although this argument has long been controversial, it has not revealed a satisfactory strategy through which the life schedule maximizes the carrying capacity.
A recent report proposed that species maximize the common objective function in both r-selection and K-selection [3]. This function provides the characteristic function of the Euler–Lotka equation—the same as in general studies [1, 4]—but it does not incorporate a parameter such as the carrying capacity. Instead, the function contains a density effect that represents the intra-specific competition from each age and state. If the density effect generates a stationary population, it indicates the carrying capacity and provides an optimal life history in K-selection. According to this model, density effects evolve various life histories not only with precocity and prolificacy but also by maximizing the population growth. This phenomenon has been observed in another study [5].
An adaptive condition for species requires not only the maximization of the population growth rate but also an evolutionarily stable strategy (ESS): no mutants can invade the population or the genetic pool. In previous research [4], the carrying capacity was considered a constraint because the objective function was assumed to monotonically decrease in terms of the density effect.
Those studies unified the population growth rate and the basic reproduction number via a characteristic function. The former is not thought to be appropriate for fitness under a saturated population. Conversely, the latter does not always become a larger population than the species maximizing the former because it does not consider the generation time. Maximizing the characteristic function maximizes the population growth rate in r–selection and maximizes the basic reproduction number in K–selection. In other words, these parameters as fitness are a result of maximizing the characteristic function and not a direct indicator of evolution. The applicability of this framework in a variable environment remains to be determined.
Evolution in a variable environment was established via sensitivity analysis [6] and Tuljapurkar’s approximation [7]. These methods have been systematized for general transition-matrix models. Recent studies focused on the effects of these structures on population dynamics in variable environments [8–11]. Each life history changes with age and has individual differences; however, it is not always reasonable to observe the growth of the physiological state with age in the field research because it is difficult to divide heterogeneity into each age structure in many cases. Therefore, biologists often apply state-structure models without age to their analysis. Essentially, researchers of empirical studies now need to consider each age structure. Evolution cannot ignore age because natural selection is believed to act with individual life histories. Multistate structured models involving age are an increasingly important area of demography and ecology [12–14].
In this study, we construct a method that addresses the adaptive life schedule in the absence and presence of a variable environment based on a continuous multi-state age-structured population model. Our method follows the general theorem for r/K-selection established by Oizumi et al. [3] and derives a more generalized control equation for the adaptive life strategy from it in a constant environment. Further, we construct a perturbation method that corresponds to Tuljapurler’s approximation in continuous models. We adopt this method for the adaptive control of heterogeneity for an optimal foraging problem in a variable environment as an analyzable example. Next, by comparing adaptive strategies in the presence and absence of a variable environment, we suggest that there exists an adaptive threshold for the variance of heterogeneity under environmental stochasticity. This study systematizes adaptive life histories in the presence of individual heterogeneity, density effects, and environmental stochasticity using the aforementioned objective function.
Our results reveal that fitness is closely related to the reproductive value. We show that characteristic functions play an important role in population dynamics even in constant and variable environments. Our model shows that heterogeneity is more likely to evolve in a variable environment than in a constant environment. Our framework will help us find a universal definition of fitness.
We developed a general model theory for biomathematics. We define the state-growth model for each trait. Suppose that y∈A⊆Rd are d-dimensional trait features characterizing each individual where A is the domain of y. The growth of each trait from age a0 to a is assumed to be described by a d-dimensional Ito-type diffusion process
Bτℓ represents the ℓ-th element of the N-dimensional Brownian motion and σjℓ(⋅) comprises
Further, gj(⋅) and Sjj′(⋅) represent the mean and covariance of j-th state growth rates, respectively.
This SDE can be interpreted as a rule for each state growth of individuals. The heterogeneity of individuals generated by the SDE is referred to as internal stochasticity to distinguish it from environmental stochasticity, which is external stochasticity [15].
For the boundary value x ∈ ∂A, each state transition rate and fluctuation term are assumed to be zero (Dirichlet condition). The age-specific fertility rate in state y is given by F(a, y) ≥ 0, and the force of mortality is assumed to satisfy
in each state because α denotes the maximum attainable age.
Let the population vector Pt(a, y), in which each individual follows the ingredients Eq (1), F(a, y), and Eq (2), be a cohort density at age a at a state y in time t. Then, we obtain the basic partial differential equation as
where the linear operator H(a, y) is given by [16]
Eq (3) implies that the cohort transitions dynamically for age a and state y at time t.
In addition, we assume that the boundary condition representing the birth law is given by
where ν(·)∈L+1(A) represents the state distribution of the neonatal population satisfying
Let pt(a) ≔ Pt(a, ⋅) represent the age-density function at time t considering a value in the trait space E = L^1^(A); further, let X ≔ L^1^(0, α;E) be the state space of the age-density functions. Then, the basic system (Eqs (3) and (4)) is formulated as an abstract McKendrick equation given as
where H(a) is a linear operator acting on E given by (H(a)f)(y)=H(a,y)f(y) for f ∈ E, and F(a) is a one-dimensional linear operator from E, given by
Suppose that the operator -H(a) generates an evolutionary system U(a, s), a ≥ s ≥ 0, on E. Then, for φ∈D(H(a)), it holds that
and the solution p is expressed as
Let β(t) ≔ pt(0) ∈ E be the density of the newborns at time t.
Substituting Eq (8) into the boundary condition of Eq (5), we have
where
Then, Ψ(a) is a one-dimensional positive operator on E, whose range is spanned by ν; the next generation operator is K=∫0αdaΨ(a). Thus, the spectral radius is given by
which is the reproduction number R0 of our system.
Let Ψ^(λ)≔∫0αdaexp{-λa}Ψ(a) and r∈C. Then, there exists a unique real root r0 satisfying the characteristic equation Λ(Ψ^(λ))=1, i.e.,
It follows from the well-known renewal theorem [17, 18] that there exist numbers C0 > 0 and η > 0 such that
where r0 is known as the dominant characteristic
and rk (k = 0, 1, 2, ⋯) are the characteristic roots of (11)
The long-term logarithmic growth rate (LLGR) of the population denoted by r¯ is defined as
where L^1^-norm ‖⋅‖X is defined as
where |⋅|E denotes the L^1^-norm of the trait space E. From the renewal theorem (9), we have r¯=r0 in a constant environment.
Let
be a linear operator on X with domain
Then, (5) can be viewed as an ordinary differential equation on the Banach space X.
where pt = pt(⋅) is a population vector taking a value in X.
Then, H becomes an infinitesimal generator of the C0-semigroup T(t), t ≥ 0, on X, and H has eigenfunctions wk as
Consider an adjoint operator H* and its eigenfunction of wk*. Let us introduce the duality pairing 〈v, w〉X between v ∈ X* and w ∈ X as
From the relationship 〈H*v,w〉=〈v,Hw〉, we have
where the domain is given by
and H*(a) is a linear operator on E* given by
The adjoint operator -H*(a) is the generator for the adjoint evolutionary system U*(a, s) = U(s, a)*, s ≥ a. It follows from (7) that
It is reasonable to define the adjoint eigenfunction corresponding to the dominant eigenvalue r0 as the reproductive value. From the adjoint eigenvalue problem Hvk=rkvk, we have the adjoint eigenvector associated with the eigenvalue r~k*~ as
where vk(0) is an arbitrary value in E.
From a stochastic perspective, transition operators U and U* are represented by a fundamental solution K(s, x → a, y) satisfying
or
(cf. [19]). Therefore, Eqs (15) and (18) can be rewritten as
Accordingly, characteristic Eq (11) becomes
This fundamental solution K(s, x → a, y) implies the transition probability of the state growth from an initial state x at age s to a final state y at age a; this is generated by Eq (1).
Using eigenfunctions, we can obtain an asymptotic expansion of the population semigroup.
where ϵ is a small positive number [20].
Further, it is easy to see that the total reproductive value V(t) ≔ ⟨v0, T(t)φ〉 satisfies
from which we have
This derivation via functional analysis is technically convenient for defining the semigroup operator using eigenfunctions; further, a stochastic interpretation of those eigenfunctions is reasonable to connect the population dynamics with the life histories of individuals. The latter interpretation is required to derive the Hamilton–Jacobi–Bellman equation involved in the adaptive control of life history, and we address this later.
To the best of our knowledge, the study of adaptive life histories using structured population models began with [1, 4]. These studies verified that maximizing the characteristic function (Eq (11)) is equivalent to maximizing the dominant characteristic root r0. Further, recent studies have extended this theorem to address internal stochasticity and density effects by adopting the stochastic control theory [3, 16].
Let us consider the general population dynamics containing the control parameter u∈U⊂Rd′, where u represents a value in the given Borel set U to control each state Xa [21]. If the d′-dimensional density effect is given by Γt∈R+d″, the general population dynamics are
Moreover, the renewal process of this system is given by
Then, if γℓ′ = γℓ′(a, y) is a weight function for each age and state, the vector of d′-dimensional density effect Γ~t~ is given by
For simplicity, H(a, y, u, Γ) is assumed to be an adjoint Fokker–Planck Hamiltonian parameterized by constant vectors u and Γ
Suppose that fertility depends on states y, u, and Γ such that
These assumptions assume that the density effects are approximated to zero or are constant in sufficiently small or stationary populations.
Here, ϕ[u] indicates that ϕ is a functional with respect to u. If u˜(a,Xa)∈U is the adaptive control of the life schedules, it should satisfy the following theorem.
Theorem 0.1Letu˜=u˜(a,Xa)be
Let**ψr [u] (Γ) be given by
and ψr0[u]u=1. We define r˜0(Γ) as r˜0≔r0u˜ as Then, we have
This theorem is easily verified because of the monotonicity of ψr [u] (Γ) with respect to r. The theorem implies that a control that maximizes ψr [u] (Γ) is equivalent to maximizing the dominant characteristic root r0(Γ) as a function of Γ (cf. [3]).
This theorem leads to two types of Let the maximized ψr [u] (Γ) be given by
One argument is related to the r selection theory that maximizes the dominant characteristic root when we choose the condition
Because Γ represents the strength of the density effects, Γ = o indicates the adaptive strategy that will satisfy the selection of r. The second argument represents the conditions in K
because the adaptive strategy in a stationary population is believed to be uninvaded by any strategy. ψ0 [u] (Γ) is essentially the basic reproductive number, and, therefore,
is necessary and sufficient for the adaptive strategy in K selection (K strategy). Γ˜ must satisfy several additional conditions, such as existence, uniqueness, and stability. The details of these additional conditions can be determined in Text A in S1 File. Although the r strategy cannot serve to conserve the exponential growth of the population in nature, it is believed to be the case that the r strategy matches the K strategy. In this case, the r strategy comprising precocity and prolificacy becomes a candidate for the adaptive strategy even in a stationary population. For example, there is a mathematical model in which intraspecific competition does not influence the control of foraging resources [3]. If ν(y) = δ^d^(x − y), our method unifies the r/K strategies via the characteristic function in Eq (27), which is matched with the consequence in the references mentioned previously.
Γ is adjusted to
assuming that each element is positive for all ℓ′:
Then, a population density P^†^(a, y) generating Γ^†^ exists and satisfies
Therefore, Γ^†^ can provide a saturated population under nonlinear population dynamics.
Let us consider the maximized function
By applying the stochastic interpretation to Eq (31), Eq (30) can be rewritten as the statistics of a diffusion process Xτ=(Xτj)1≤j≤d as
where Ey[·] denotes the expectation of the probability measure of Xτ in Xa = y. This representation is called the Feynman–Kac formula, and it is well known in stochastic analysis [19]. Eq (32) is called the value function in the control theory [21]. The diffusion process Xτ=(Xτj)1≤j≤d satisfies the following stochastic differential equation (SDE):
The SDE is given by Eq (1) parameterized by u and Γ, and it can describe the growth process of each state from age a to u in both trivial (Γ = 0) and nontrivial (Γ = Γ^†^) equilibrium points.
Thus, vr [u] (a, y, Γ), the solution of the Dirichlet problem
provides a statistical representation of the corresponding diffusion process called the Feynman–Kac formula [19, 22]. The adjoint Hamiltonian is given by
The stochastic interpretation is appropriate for describing the adaptive life history and corresponding population dynamics for the following two reasons. (1) To reveal that the fittest dynamics are generated by the optimally controlled life history of individuals. (2) To derive the main equation in this study from the central principle of optimality efficiently.
According to the optimal control theory, adaptive strategies must follow a basic property called Bellman’s principle (or the principle of optimality):
“an optimal strategy has the property that whatever the initial state and initial control are, the remaining control must constitute an optimal strategy with regard to the state resulting from the first strategy” [23].
The following relationship is derived based on this
where 0 ≤ a0 ≤ a ≤ α. This relationship implies that the adaptive control from a0 to a in the terminal condition v˜r(a,y,Γ) is consistent with the control of this function from a0 to α, and it leads to
(see Text B in S1 File). This equation is significant in control theory and is called the Hamilton–Jacobi–Bellman (HJB) equation. From the basic theorem of the adaptive life schedule, the adaptive strategy u˜r(a,y)|y=Xa in r selection (r strategy) is obtained using Eq (28), and the K strategy u˜K(a,y)|y=Xa is given by Eq (29).
Based on Eqs (28) and (29), Eq (34) is simplified as
Thus, we obtain an equation for which the adaptive strategy is satisfied in a constant environment. Eqs (28), (29) and (36) contain and are more general than the result of [3] because they account for reproductive controls. Moreover, these equations reveal that adaptive control depends on the state distribution of the neonatal population ν(y) via r˜0 or Γ˜. Accordingly, the equation above indicates that individual life histories evolve to maximize the reproductive value function (Eq (32) at age zero) in a constant environment.
The previous sections revealed a parameter that maximizes the adaptive life history in a constant environment. This section presents the population dynamics behavior under a simple stochastic environment.
Although there are several assumptions and candidates for statistical noise as external stochasticity, we simplify environmental stochasticity as white noise parameterized by a and y Wt(a,y).
for all t > 0. Bt(a,y) denotes the Brownian motion parameterized by a and y. Consider that a population vector under external stochasticity Ptε(a,y) follows the stochastic partial differential equation
where ε denotes a sufficiently small positive constant that represents the strength of external stochasticity. Because it is difficult to compute a strict value of an LLGR involving external stochasticity, we apply a perturbation method to ε to calculate its approximate value, such that
We introduced the derivation of the second-order approximation of LLGR in Eq (40). The population Hamiltonian vector, Hamiltonian, and noise functions are simplified to avoid computational complexity as
Let us consider the following variation of the constants
The semigroup T(t) is defined by Eq (21). With Eq (41) and Ito’s formula for the multiple stochastic integral [24], a perturbation of the population vector is found by computing iteratively.
Introducing a new operation symbol
each Qm (t) φ is deformed as
where hm(x) denotes the Hermite polynomial
The last row in Eq (43) is derived from the following formula [24]:
The arbitrary constant of the adjoint eigenfunction is set as
If the population vectors in the presence and absence of external stochasticity are close to each other, Ptε≈Pt0=Pt(ε⪡1), the perturbation expressed in Eq (43) provides an accurate approximation. With this assumption, a ε-specific mean LLGR r¯E(ε) is represented by substituting Eq (43) into Eq (23) such that
For simplicity, suppose that the initial population is the eigenfunction corresponding to the 0–zeroth characteristic root.
By expanding Eq (44) into a Taylor series at ε = 0, the growth rate becomes
Let us consider the mean growth rate in an environment comprising sufficiently small disturbances such that the third- (or higher-) order terms in ε can be truncated. The second term on the right-hand side is zero in the mean growth rate because of the statistical property of the fluctuation term (cf. Eq (38)). Accordingly, the key point is the estimation of the second-order term in Eq (45). One of the pieces composing the second-order term is computed as
Hermite polynomials are orthogonal with respect to Gaussian measure
i.e., the term becomes statistically zero.
Similarly, the other component of the third term in Eq (45) is computed as
After combining the components (Eqs (47) and (48)), the second-order approximation of the LLGR becomes
This approximation is similar to the Tuljapurkar approximation [7]; however, it differs in several aspects. For instance, the deviation term corresponding to the original Tuljapurkar approximation is described by a sensitivity matrix. In this continuous version, statistics concerning the diffusion process Xa account for the term. One important point is that the second term on the right-hand side of the equation above incorporates eigenfunctions. As described previously, the adjoint eigenfunction serves as an objective function to determine the adaptive strategy. This characteristic suggests that an adaptive species in a variable environment does not always maximize identical functions in a constant environment. That is, we may find another adaptive strategy u* as
where the arbitrary constant is set to
The previous section revealed that the effect of external stochasticity on population growth is represented by the eigenfunctions corresponding to the dominant characteristic root in the mean environment. We use a specific mathematical model that is analytically solvable to examine the contribution of internal stochasticity to external stochasticity.
Let us consider the role of internal stochasticity in external stochasticity. We construct a mathematical model that compares the LLGRs on a group of inhomogeneous growth rates with those of a group of homogeneous growth rates in a variable environment (cf. Fig 1). This figure illustrates the concept of a simple model. This model verifies whether the variance in size growth σ1 increases with the LLGR r¯E(σ1) for positive values of ε.
Fig 1 Variation in individual size growth in a variable environment.
The model aims to estimate the existence of the adaptive control of internal stochasticity against external stochasticity. As indicated in the aforementioned analyses of matrix models based on empirical data, if organisms control their growth rate statistics, there exists an adaptive strength of heterogeneity.
When Xa∈R+ is the size at age a ∈ [0, ∞), as an effect of internal stochasticity, we assume that the heterogeneity of the individual size growth rate is
Heterogeneity is generated by the fluctuation in the second term on the right-hand side of Eq (50). The SDE provides a geometric Brownian motion that grows exponentially with the fluctuation. Suppose that mortality is constant. Then,
Fertility is assumed to be an allometric function in size
This life history generates the following Hamiltonian and adjoint Hamiltonian.
respectively. Assuming all neonates have identical state x
eigenfunction wr satisfies
Substituting an ansatz
into Eq (55), the equation is converted into a Fokker–Planck equation
which gives the probability density function of the geometric Brownian motion in Eq (50). The probability density function is then given by the logarithmic normal distribution
Therefore, the eigenfunction is
Because the size growth rate follows an age-homogeneous Markovian process (Eq (50)), this adjoint function does not depend on age.
Therefore, we have
Then, the adjoint eigenfunction follows the adjoint equation
This equation is explicitly solvable using the following
which gives
Because the function above can compose the characteristic equation
the dominant characteristic root is computed as
By the definition of 0 < ρ < 1, the characteristic root indicates that internal stochasticity has a negative effect on population growth in a constant environment
That is, it is nonadaptive for species to have heterogeneity under the condition 0 < ρ < 1.
Substituting the dominant characteristic root Eq (62) into Eq (60), the functional becomes
Hence, the arbitrary constant is determined to be
and the adjoint eigenfunction corresponding to the dominant characteristic root is
In this case, the adjoint eigenfunction corresponding to the dominant root matches the fertility function.
Combining the previous steps, the deviation term is given by
This deviation term diverges to infinity in the absence of internal stochasticity
Thus, it is reasonable to consider the effect of higher orders of ε on this divergence. However, this consequence suggests the significance of heterogeneity in the persistence of species in minimally variable environments. This property contrasts with the effect of internal stochasticity on the dominant characteristic root (Eq (62)). Substituting Eqs (62) and (64) into Eq (49), the LLGR approximates
The LLGR represents a monotonically increasing function with respect to the mean size growth rate b1,
This point may appear to be trivial, yet it is notable that the deviation term monotonically decreases in b1. Further, rapid growth may reduce the risks inherent to variable environments.
The hHeterogeneity of the size growth rate reduces the mean dominant characteristic root and causes the risk of extinction from the variable environment. By computing Eq (65) in terms of σ and ε, we find the adaptive heterogeneity of the size growth rate for each ε (see Fig 2). This figure shows the existence of an adaptive value in σ1. Each ε representing the strength of external stochasticity has a unique adaptive value of σ1 that maximizes the dominant characteristic root r0; ε increases with an adaptive value σ1. This result suggests that species require greater heterogeneity in more variable environments. The parameters are b1 = 0.6, x = 0.01, μ0 = 0.1, f0 = 1.0, and ρ = 0.4.
Fig 2 Adaptive heterogeneity under two-fold stochasticity.
Fig 2 illustrates that adaptive heterogeneity increases with environmental variability. The numerical analysis suggests that species evolve to yield heterogeneity in variable environments. This viewpoint corroborates conventional interpretations of the necessity for biodiversity.
Based on Eqs (50)–(52), we consider a species utilizing different resources (R1 and R2). The specialist utilizing R1 uses the size growth rate in Eq (50) and that of a specialist utilizing R2 is
Ba1 and Ba2 are independent Brownian motions. Then, we assume that b1∈R+≥b2∈R (b2 could be negative), σ1 > σ2 ≥ 0, that is, choosing R1 implies a higher risk and growth rate expectation than choosing R2. Conversely, choosing R2 under the same conditions confers another risk—that individuals have lower survival until they reach maturity than when choosing R1 because of the slower average growth rate. Therefore, individuals should find their adaptive risk by hedging u˜(a,Xa)∈[0,1] in accordance with each population size under the following growth rate (cf. Fig 3).
Fig 3 Adaptive resource utilization in external stochasticity conditions.This model is referred to as the two-resource utilization model and it generates the following adjoint Fokker–Planck Hamiltonian.
Fig 3 illustrates the concept of the adaptive resource utilization model. A resource R1 provides the high size growth rate b1 on average; however, the risk σ1 is also high. Conversely, R2 is low risk σ1 > σ2 and has a low size growth rate on average, i.e., b1 > b2. The species maximizes its LLGR by optimizing the utilization of both resources. Then, we verify that the existence of external stochasticity evolves different adaptive utilizations from that in a constant environment.
In this model, finding the adaptive utilization is analogous to generatinge the optimal size growth curve with heterogeneity. This growth curve maximizes Eq (57) following our framework. Consequently, individuals adopting the adaptive allocation strategy compose the fittest species by maximizing the LLGR under r-selection.
Because the reproductive value is independent of age in this model, the value function (Eq (32)) also does not depend on age, such that
From the Bellman’s principle Eq (33), the value function has the following decomposition for arbitrary age a:
The equation above is deformed as
Using the same process as that used for the derivation of the general HJB equation (see S.2.), adopt the Feynman–Kac formula [19, 21] into the equation
Take the limit as a tends to zero such that
Then, we have
The equation above implies that the adaptive control should provide an extreme
for all x. This necessary condition leads to the following relationship between adaptive utilization and the adjoint function.
Thus, the control is independent of age, which is called stationary control in control theory. Substituting the adaptive control condition into the adjoint Hamiltonian
we can derive the adjoint eigenfunction of the adaptive life history from the same ansatz, as in Eq (60).
From Eq (72) and the function above, adaptive utilization is computed as
which is identical to the strategy in [16], and it is known as constant value control. It indicates that R2-specific utilizers do not evolve. Because adaptive utilization is constant in constant environments, finding another utilization constant u* that maximizes the LLGR in a variable environment implies that another adaptive utilization exists, even if the constant is not optimal control. Suppose that utilization is always constant, and that v* becomes the adaptive strategy for twofold stochasticity. The utilization constant specific LLGR r¯u considers whether the variable environment selects a life history that favors heterogeneity as adaptive. Because the utilization rate does not depend on age or size, we can consider a specific LLGR with the following change of coefficients in Eq (65).
Solving r¯u numerically, we find that a variable environment favors heterogeneity, as suggested in the previous section (Fig 4). This figure illustrates adaptive resource utilization with respect to ρ for several values of ε. Although ρ represents the scaling exponent of fertility that denotes a measure of risk aversion in a deterministic environment, the adaptive utilization of risk appetite in the presence of external stochasticity exists in the domains of greater and smaller values of ρ. The domain of the adaptive strategy utilizing both resources narrows as ε increases. This consequence is linked to the relationship between the adaptive value of σ1 and ε in Fig 1. The parameters are b2 = 0.5, σ1 = 0.8, σ2 = 0.005, and ε = {0, 0.1, 0.3, 0.6, 0.9}; the others are the same as in Fig 1.
Fig 4 Adaptive resource utilization under external stochasticity.
The scaling exponent ρ represents a risk appetite index in economics; small values favor risk aversion. Adaptive strategy Eq (73) represents an identical interpretation of the exponent to that in economics. However, under external stochasticity ε ≠ 0, minimal internal stochasticity does not become adaptive for small values of ρ.
This study attempted to construct a systematization of the optimal life schedule problem and population dynamics using the eigenfunction expansion of a structured population model. Our perturbation method was inspired by Tuljapurkar’s approximation; however, our model is based on mathematical models of life scheduling that contain internal stochasticity (e.g., SDE). This change provides a theoretical basis for the argument that species pose environmental variability. By applying the framework in this study to this argument, we found that the optimal parameters reduce the risk of external stochasticity and increase the LLGR. Further, the conventional adaptive life history in constant environments can be found using the HJB equation derived from the adjoint equation. If we regard the ESS in stable populations as an adaptive strategy, then using the HJB equation with additional parameters to represent the magnitude of the density effects on state in a stationary population can provide an adaptive life history under intraspecific competition.
The framework of this study helps reveal what evolution maximizes. In a constant environment, this framework extends the consequence in [3]: both adaptive strategies in the presence and absence of density effects maximize a common function. Further, although adaptive strategies under variable environments are less simple than those in constant environments are, this study shows that the effect of external stochasticity is closely related to eigenfunctions and Tuljapurkar’s approximation. The second-order perturbation used in this study yielded trade-offs between the mean dominant characteristic root and the corresponding eigenfunctions via the LLGR. Given this relationship, the effect of internal stochasticity on the population growth may differ from its effect in a constant environment. Thus, as shown in the analysis of the specific model, the same adaptive strategies are not always used. External stochasticity needs to be treated as a different type of selection pressure than internal stochasticity and density effects. Therefore, deterministic models approximated by the averaged environment often overlook the essential adaptive strategies.
The specific model showed that deviation in the size growth rate buffers the reduction in LLGR caused by the variable environment. There is a trade-off between the decrease in the mean characteristic root and the buffering effect of a variable environment with internal stochasticity. This determines the adaptive heterogeneity of size growth, where the length is proportional to the magnitude of external stochasticity. These consequences support the premise that adaptive utilization prefers high-risk resources in a resource-utilization model in a variable environment. Compared with a constant environment, the domain of the allometric exponent wherein species evolve with risk aversion under constant environments changes to a risk-taking strategy. Despite the small exponent value of a species, which indicates the brittleness of internal stochasticity in the deterministic LLGR, risky resources in a variable environment are selected. Therefore, these specific models appear to provide a theoretical basis for the conventional argument that individual heterogeneity is necessary for living in a variable environment. The mean rapid size growth in risky resources can be interpreted as having an advantage in terms of a small exponent value because fast growth statistically reduces the risk of external stochasticity. Considering that precocious species, such as mice, have short lifespans, this interpretation may be related to the short lifetime of organisms in variable environments [10, 11, 25]. However, in such a simple model, this interpretation requires careful consideration because the deviation term of the LLGR does not depend on mortality. Despite this simplification, our framework links empirical studies of evolution that pertain to life histories to various theoretical studies of structured population models.
The perturbation method in this study also avoids the mathematical complication of external stochasticity at the expense of biological correctness; incorporating these features remains an open problem. For instance, all cohorts should monotonically decrease with age; however, setting white noise as external stochasticity violates this rule. There is a limit to this study. Eq (40) can be interpreted as the fluctuation of mortality from external stochasticity; however, white noise neither correlates with each age-containing parameter nor ensures the positivity of mortality. This assumption is only for the sake of mathematical simplicity because it allows us to assume that external stochasticity alters the mean state growth rate and the fluctuation term from internal stochasticity or both. In this case, note the treatment of the derivatives in the noise function; these assumptions can obey the aforementioned biological rules because of the conservation law in the continuity equation. These noise functions are thought to complicate the problem and require considerable mathematical discussion. In addition, candidate stochastic processes are believed to vary such that the SDE can be defined by Ito’s integral, Stratonovich’s integral, and others [26, 27]. If we choose a noise function that does not have Markovian properties, the approximation of LLGR may not correspond with our results.
Disregarding the configuration method of each stochastic process, this study was conducted under the premise that all stochastic processes are assumed to be Markovian, which is an assumption that has been accepted by many ecologists. Structured population models have various versions including age-, size-, and stage-structured models; Tuljapurkar’s approximation appears to work well in size- and stage-structured models that ignore cohort information. However, individuals and cohorts are essential elements when considering evolution in variable environments. Many empirical studies based on models that exclude cohort dynamics suggest a correlation between the transition rate and environmental variability [8]; however, they cannot clarify the strategy by which every individual’’s life history reduces the risk of external stochasticity. On the other hand, these empirical studies suggest that the vital rate, which is important for adaptive strategy, is robust against environmental changes [9]. This suggestion imposes an important requirement on theoretical studies of life history evolution in a variable environment. Theoretical studies based on cohort dynamics should also consider this requirement. As mentioned in the Introduction, a twofold stochasticity perspective of cohort dynamics is necessary to understand the effect of stochasticity on life histories. Future transition matrix models must consider the age structure to understand how organisms oppose risk in a variable environment in their life history.
Thus far, theoretical research on the evolution of life history has been focused on an individual [28]. Within a lifespan, the strategy of maximizing the basic reproduction number is considered to be adaptive. The idea was the same in a variable environment [29]. The drawback of maximizing the basic reproduction number is that the generation time is not considered. Therefore, it does not always match the maximization of r0 and LLGR. In r–selection, the maximization of basic reproductive number may not be the optimal solution. Our framework overcomes that problem.
As shown in the analysis of the specific model in this study, internal and external stochasticity yield the diversity and extinction of organisms. Despite the simplicity of the assumptions of stochasticity, this study quantitatively demonstrates that heterogeneity decreases risks associated with a variable environment. Further, this result suggests that the existence of adaptive heterogeneity maximizes the population growth. Because many organisms are believed to have various adaptive strengths of traits, the diversity on an ecological scale may also occur. Eqs (28), (29) and (36) link the evolution of life history to population growth under internal stochasticity. Eq (49) connects the life history with the effect of external stochasticity via eigenfunctions. In r selection, an adaptive strategy must optimize not only the basic reproductive number but also the generation time. An adaptive strategy in K selection must generate density effects that prevent a stationary population from being invaded by other strategies. A previous study [3] posited that adaptive strategies in both r and K selection were identical via the common HJB equation, and this provides adjoint eigenfunctions.
On the other hand, the density effect from other states will depend on the current state. In addition, in terms of fertility, parental status is generally thought to affect the initial status of the offspring. In this study, these state dependences were ignored and assumed to be constant. Eliminating these assumptions will allow us to express more realistic intraspecific competition.
Consequently, this study shows that r and K selections and external stochasticity evolve different phenotypes; these selection pressures are independent of each other. In r selection under a constant variable, our simple model shows that the heterogeneity should decrease because it decreases the expectation of the characteristic function. In the K selection, the previous study demonstrated that the evolution of heterogeneity depended on how density effects operated in life history [3]. In r selection under the variable environment, the homogeneity poses a high risk of extinction. The last one is caused by a trade-off between the mean growth rate and its variance in population dynamics. However, these results also suggest that the consequences of evolution in life history arise from optimizing a common factor, i.e., the reproductive value, in each habitat. To prove this, we must examine whether a life history adaptive strategy in more complicated environments (i.e., containing both density effects and external stochasticity) is explained by the framework developed in this study. The choice of the density effect and the definition of background noise (including non-Markovian) will generate numerous evolutions concerning heterogeneity in life history. Studies on these themes will find more sophisticated concepts of fitness. We hope that this research will be one of the cornerstones for future research.
We are deeply grateful to Yuki Sugiyama for his helpful comments and advice on this study; Nobuhiko Fujii, Kensaku Kinjo, Youichi Enatsu, Kota Hattori, Tetsuya J. Kobayashi, and Akira Sakai for their mathematical advice; and Kumiko Oizumi, Kunihiro Aoki, Hiroko Oizumi, and Ryuichi Kaneko for their support and encouragement.
All relevant data are within the manuscript and its Supporting information files.
This research was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI (Grant number 20K14368) and the Ministry of Health, Labour and Welfare (Grant number 20AA2007).