# Optimal allocation in annual plants
with density dependent fitness

Sergiy Koshkin, University of Houston-Downtown, Houston, TX, US
email<EMAIL_ADDRESS>
Zachary Zalles, Rice University, Houston, TX, US,
Michael F. Tobin, University of Houston-Downtown, Houston, TX, US
Nicolas Toumbacaris, Emory University, Atlanta, GA, US
Cameron Spiess, Keene State College, Keene, NH, US

###### Abstract

We study optimal two-sector (vegetative and reproductive) allocation models of
annual plants in temporally variable environments, that incorporate effects of
density dependent lifetime variability and juvenile mortality in a fitness
function whose expected value is maximized. Only special cases of arithmetic
and geometric mean maximizers have previously been considered in the
literature, and we also allow a wider range of production functions with
diminishing returns. The model predicts that the time of maturity is pushed to
an earlier date as the correlation between individual lifetimes increases, and
while optimal schedules are bang-bang at the extremes, the transition is
mediated by schedules where vegetative growth is mixed with reproduction for a
wide intermediate range. The mixed growth lasts longer when the production
function is less concave allowing for better leveraging of plant size when
generating seeds. Analytic estimates are obtained for the power means that
interpolate between arithmetic and geometric mean and correspond to partially
correlated lifetime distributions.

Keywords: allocation schedule, fitness, evolutionary stable strategy, eco-
evolutionary feedback, density dependence, temporal variability, recruitment
survival, seed yield, bet hedging, random lifetime, diminishing returns, mixed
growth, graded allocation, control-affine system, singular control

MSC: 37N25 49K15 49K30 92C80 49N90 92D15

##  1 Introduction

Models of allocation to growth and reproduction in plants based on the optimal
control theory were pioneered by Cohen [6] (in the continuous case, by Denholm
[12]), and developed in a variety of ways by many authors [10, 28, 34, 44, 46,
49]. Approaches and results are summarized in the survey paper [36], and the
book [39]. These models divide plants into two sectors, the vegetative
(leaves, stems, roots) and the reproductive (flowers, fruits), and the model
optimizes the allocation of net production between the two. Namely, it
prescribes how to vary the fractions of net production allocated to the
vegetative and to the reproductive growth over plant’s lifetime to maximize a
fitness measure. Models with more than two sectors, that prescribe shoot and
root allocation separately [25], or add storage [34], have been studied as
well.

There are two sorts of assumptions about plants’ lifetimes in allocation
models of this type. In the more common deterministic models the lifetime is
assumed to be fixed and the reproductive mass (assumed to be proportional to
the seed yield) at the end of it is taken as the fitness measure. Two sector
models then typically predict as optimal bang-bang behavior, which is
allocation to only vegetative growth followed by allocation to only
reproductive growth, with the switching time determined by the length of
lifetime and the production rate [12, 46]. Models with a random distribution
of lifetimes are more realistic but rare [1, 10, 13, 28, 31], and typically
assume that either geometric or arithmetic mean of seed yield is maximized.
For geometric mean maximizers, a novel phenomenon is that the optimal behavior
may no longer be bang-bang, but rather there is also a period of mixed growth
when vegetative and reproductive organs grow simultaneously.

The main goal of this paper is to investigate effects of density dependence
and temporal variation on the choice of fitness measure and the resulting
optimal allocation schedules for annual plants. Density dependence is a simple
example of eco-evolutionary feedback, where selection patterns in the
population are affected by growing density of the population itself [38].
There is an ongoing discussion on what, if anything, is maximized by natural
selection, especially when random fluctuations are taken into account [14, 16,
40], and appropriate fitness measures often depend on the assumed type of
feedback. Fisher’s Malthusian parameter $r$ is a good fitness measure only for
populations growing in a vacuum, with no appreciable density effects.
MacArthur and Wilson termed selection under such conditions $r$-selection, and
showed that in environments saturated with organisms $K$-selection takes place
instead [33], where $K$ is the “carrying capacity” of the environment, or,
more precisely, the stable population size. For many species $K$-selection can
be associated in deterministic models with maximizing the lifetime
reproductive success $R_{0}$ [14], which can be measured by seed yield for
annual plants. The oft-used geometric and arithmetic mean maximization in
fluctuating environments also correspond, roughly, to $r$ and $K$ selection,
respectively.

However, already Pianka suggested that in real environments the two extremes
are compromised giving rise to a continuum interpolating between $r$ and
$K$-selectors [37]. Moreover, if the environment is fluctuating then both
types of selection might be contributing to the final outcome, and the
environment itself might be affected by the selection taking place. More
recent modeling suggests that fluctuating environments and density dependent
feedbacks lead to transitions between fitness measures such as $r$ and
$R_{0}$, and to maximizing mixed or hybrid measures that interpolate between
the two [3, 4, 14, 17, 29]. Competition between different genotypes means that
only relative difference in reproduction across competing genotypes matters,
which implies diminishing returns on selective advantage from its increase
[16]. In agreement with the method of evolutionary stable strategies, such
fitness measures reflect invasibility of a trait and/or the chance of its
fixation in the long run [17, 29].

Furthermore, the notions of $r$ and $K$ selection are linked to the classical
logistic model, which has been criticized for including obscure
phenomenological parameters and unrealistic predictions at the equilibrium. A
recent model of nest site lottery, proposed in [2] and developed in [3, 4],
gives a deeper insight into actual demography and differences in selection for
unsuppressed and suppressed growth, which stand behind the intuitions for $r$
and $K$ selection, respectively. It provides an example of a mechanism of
selectively neutral growth suppression acting on the recruitment of juveniles
(seeds and seedlings in the case of annual plants), and shows that
transitional fitness measures should reflect more than just a difference in
averaging. Juvenile mortality directly reduces reproductive value of seeds,
and not necessarily proportionally. While availability of nest sites as a
growth suppression factor is not applicable to plants, analogous mechanisms
operate on seeds and seedlings as well.

As we argue below (Section 1.1), both mixed averaging and recruitment survival
effects can be incorporated into a non-linear function of reproductive mass
$L(y)$ that enters a fitness measure of the form $\mathbb{E}[L(y)])$, where
$\mathbb{E}$ is the expected value averaged over individual lifetimes. We call
$L(y)$ the fitness function. When the number of surviving seeds is
proportional to the total reproductive mass, i.e. when recruitment survival
can be neglected, one can take $L(y)=y$ for arithmetic and $L(y)=\ln y$ for
geometric mean maximization, and those are the only two cases typically
considered in the literature. But when non-linear effects of recruitment
survival are taken into account, a much broader range of $L(y)$ is realized,
from superlinear to sublogarithmic.

The role of $L(y)$ is muted in deterministic models, because maximizing $y$ is
equivalent to maximizing any monotone function of it, but its effect can be
rather dramatic in allocation models with random variation, because
$\mathbb{E}[L(y)]\neq L(\mathbb{E}[y])$. Mathematically, fitness functions are
analogous to utility functions in economic models, although their
interpretation is quite different from describing the subjective value of
returns on investment. One can think of $L(y)$ as the “evolutionary utility”
of reproductive output $y$. In particular, linear and superlinear (convex)
functions represent “gambler’s utilities” favoring high risk high reward
strategies of delaying reproduction, while concave functions lead to more
conservative behavior that Slatkin called “hedging one’s evolutionary bets”
[36, 41, 43].

In addition to general fitness functions, we also consider general production
functions assumed to be concave, instead of the typical linear ones, to
reflect diminishing returns. In plant biology, concave production functions
are natural as a consequence of Liebig’s law of limiting factors [19]. For
example, self-shading of leaf canopies limits returns on increasing the leaf
mass, other limiting factors include availability of minerals and nutrients.

###  1.1 Reproduction and fitness

The fitness measure of our model takes into account both recruitment survival
of seeds to seedling establishment, and temporal variability of lifetimes
within a season. The two factors have to be taken into account separately
because the allocation model applies only after plants begin to rely on
autotrophy rather than on maternal resources for growth (as is typical). Since
the model’s lifetime distribution only covers established plants the seed
yield as a measure of reproductive success has to be corrected for recruitment
survival. In the long run, advantages and disadvantages of allocation
schedules are sensitive not to the total number of seeds produced, but to the
number of them that will live to reproduce, in the next or subsequent seasons.
Indeed, the tradeoff between the risk of delayed reproduction and potentially
larger seed yield may or may not be worthwhile depending on how many extra
seeds will likely grow into reproducing adults. Even if we assume that the
fraction of the overhead costs (flowers, pedicels and the like) is size
independent and the total seed yield is proportional to the reproductive mass
$y$, it is not the total yield that determines reproductive success, but
rather the effective seed yield $S(y)$. The effective seed yield is the
expected number of seeds that survive seed and seedling mortality, and are
recruited into the adult population that may reproduce ($S$ stands for
‘surviving’). We model $S$ as a function of $y$ only, which is a
simplification that can only be justified in special circumstances. In
particular, it does not take into account that survivability may vary in the
course of a season.

When the effective seed yield is averaged over individual lifetimes to produce
a fitness measure, an additional non-linearity may be introduced. This is
familiar from the case of geometric mean, where the averaging transformation
is $A_{0}(z)=\ln z$. More generally, mixed fitness measures correspond to
power means with $A_{\alpha}(z)=z^{\alpha}$, where $0<\alpha\leq 1$ is
determined by the correlation between individual lifetimes in the population,
as we show below. The overall fitness function is then the composition of the
two non-linearities, i.e. $L(y)=A_{\alpha}(S(y))$, perhaps scaled and/or
shifted by a constant for convenience. Thus, both recruitment survival and
averaging contribute to the non-linear dependence of the fitness measure
$\mathbb{E}[L(y)]$ on the reproductive mass.

####  1.1.1 Averaging

We start with the averaging as the simpler of the two factors. The
conventional wisdom of geometric mean maximization goes back to a simple two
allele haploid model of Dempster [11], which showed that the allele with the
larger geometric mean gets fixed with greater probability. In a number of
papers [1, 5, 6], Cohen and co-authors reached the same conclusion and
explored its implications for the behavior of annual plants. More complex
models confirmed the principle, but only under the assumption that the
correlation of lifetimes between individuals of the same allelic type is near
perfect (and cross-correlation among types is near zero) [16], e.g. because
most of them die near the end of the season. Lewontin-Cohen [31] and later
Gillespie [18] also considered the other extreme, when the correlation is near
zero, and found that the appropriate fitness measure is then the arithmetic
rather than the geometric mean of seed yields.

Quantitatively, the effect of temporal variability on the fitness measure was
studied by Frank and Slatkin [16, 17], and more recently in [29]. In general,
the genotype with the greater probability of trait fixation has the larger
value of $m-\rho\frac{\sigma^{2}}{2m}$, where $m$ is the mean, $\sigma^{2}$
the variance, and $\rho$ the correlation coefficient between lifetimes. For
$\rho=0$ we get the arithmetic mean, and for $\rho=1$ a mean-variance
approximation of the geometric mean. As $\rho$ varies between $0$ and $1$ the
expression $m-\rho\frac{\sigma^{2}}{2m}$ can be seen as the mean-variance
approximation of power means
$M_{\alpha}:=(\mathbb{E}[y^{\alpha}])^{\frac{1}{\alpha}}$
($M_{0}:=\exp(\mathbb{E}[\ln y])$, the geometric mean), where $\mathbb{E}$ is
the mathematical expectation. Indeed, $M_{\alpha}\approx
m-(1-\alpha)\frac{\sigma^{2}}{2m}$, so $\alpha=1-\rho$. Thus, we can expect
power mean maximization when plant lifetimes are less than perfectly
correlated.

If we assume that $S(y)$ is proportional to $y$, i.e. the same proportion of
seeds survives regardless of their total number, one can take
$L(y)=A_{\alpha}(y)$ because a constant coefficient of proportionality does
not affect optimization. The corresponding fitness measure, which serves as
the payoff functional in optimization, will then be of the form
$\mathbb{E}[A_{\alpha}(y)]$, where $A_{\alpha}(y):=y^{\alpha}$ for
$0<\alpha<1$ rather than $A_{1}(y):=y$ or $A_{0}(y):=\ln y$.

According to Amir and Cohen [1, 6], “linear” (i.e. arithmetic mean) maximizers
are favored in large populations and stable environments, while “logarithmic”
(i.e. geometric mean) maximizers are favored in smaller populations and
harsher environments. It is not hard to see the parallel between logarithmic
and linear maximizers, and Pianka’s $r$ and $K$-selectors, respectively.
Transitional allocation patterns produced by environments that fall between
this dichotomy have not, to our knowledge, been studied in the plant science
literature. Biological or ecological processes reducing the correlation among
lifetimes for a cohort of seeds will shift the selective advantage for
allocation from geometric mean maximizers toward arithmetic mean maximizers,
but not all the way. For example, some plants produce cohorts of seeds which
germinate over multiple seasons due to delayed germination of a subset of the
seeds [45]. Seed dormancy for more than a year (seed bank) converts seasonal
variability into cohort variability. Similarly, correlation among plant
lifetimes could be reduced by spatial variability in resources due to the
differential impact of resource availability on survival of spatially
dispersed plants [8]. Disturbances like fire could increase or decrease the
correlation depending on the spatial variability of fire induced mortality
[35], thus shifting selection toward geometric or arithmetic maximizers,
respectively.

####  1.1.2 Recruitment survival

As we saw, different averaging leads to a very narrow class of fitness
functions. But it broadens considerably when non-linear effects of seed
survival are also taken into account, i.e. when the averaging is applied not
to $y$ directly but to $S(y)$.

Ecological processes resulting in nonlinearity between seed yield and
effective seed yield are very heterogeneous and complex. Among other things,
seed survival is affected by seed dormancy, size, dispersal, physical or
chemical defenses, and density dependent inhibition by other seeds, seedlings
or adults, with these mechanisms selected for or against in concert with the
allocation schedule. Determination of $S(y)$ is further complicated by the
fact that seeds germinate in different years due to seed banks, and at
different sites due to dispersal. Nonetheless, field studies of seed
germination suggest [23, 32] that seed and seedling density has some direct
effects on germination, both positive and negative, that can be captured by
$S(y)$. If the proportion of germinating seeds that successfully establish
themselves decreases as seedling density increases we would expect a concave
dependence of effective seed yield on total seed yield.

More generally, effects of selectively neutral juvenile mortality can be
captured by $S(y)$. Even simple known mechanisms that lead to mixed fitness
measures, such as the nest site lottery [2, 3] and its generalization, the
host lottery that models parasitism [4], have close analogs for annual plants.
For example, a limited number of sites for establishment [9] would result in a
fixed carrying capacity of the environment (in the mechanistic sense of [21]
rather than the phenomenological sense of the logistic growth model). The
resulting function $S(y)$ will be linear for small $y$, but will saturate at
the carrying capacity. A more flexible mechanism of the same nature is
germination suppression to reduce competition when the seed or seedling
density is high [23, 32], which leads to more gradually saturating curves.
Intraspecific competition among seedlings for resources, or herbivores and
pathogens attracted to greater densities of host plants are also potential
mechanisms producing a concave dependence [20].

Alternatively, a convex dependence can result from predator satiation.
Predator satiation has been invoked to explain masting in perennial species
[27], but modeling suggests that this type of selective pressure could favor
temporal clumping of seed production in annuals as well [7]. Granivores often
consume a lower proportion of seeds as the amount of seeds available
increases, leading to the inverse density dependence and thus a convex shape
of $S(y)$ [22]. Similarly, synchronized flowering or fruiting may improve
efficiency of pollination and seed dispersal [27]. For example, larger seed
crops were associated with more efficient dispersal by rodents in a forest
system [48].

As we saw, when $S(y)$ is linear one can take $L(y)=y^{\alpha}$ for
$0<\alpha\leq 1$ depending on the type of averaging alone. On the other hand,
one can expect initially convex dependence that saturates at larger values of
$y$, such as given by functions of the form $S(y)=S_{\max}e^{-Cy^{\alpha}}$
with $\alpha<0$, from a combination of parasite/disperser suppression with a
carrying capacity limit, see Figure 1(a). Composed with the geometric mean
averaging $A_{0}(y)=\ln y$ this gives $L(y)=-Cy^{\alpha}+\ln S_{\max}$, where
the multiplicative and the additive constants can be dropped for the purposes
of optimization, leaving just $L(y)=-y^{\alpha}$.

(a) (b)

Figure 1:  Possible graphs of effective seed yield as a function of
reproductive mass: a) $S(y)=S_{\max}e^{-Cy^{\alpha}}$ with $\alpha<0$; b)
$S(y)=S_{1}y^{\gamma}$ with $\gamma>1$.

One can also consider $S(y)$ that remain convex over a feasible range of $y$
values, such as $S(y)=S_{1}y^{\gamma}$ with $\gamma>1$, for a suitable
normalization constant $S_{1}$, see Figure 1(b). Composed with the power
averaging $A_{\alpha/\gamma}$ they produce the fitness function
$L(y)=S_{1}^{\alpha/\gamma}y^{\alpha}$, where the constant can again be
dropped. It is convenient to unify the $\pm y^{\alpha}$ cases by taking
$L(y)=\frac{y^{\alpha}}{\alpha}$ for all real $\alpha\neq 0$, which
automatically takes care of the sign and simplifies some formulas. We refer to
$\alpha$ as the fitness index. Of course, much more general functions $S(y)$,
and hence $L(y)$, are possible, and some of our general results apply to them.

###  1.2 Model description and main results

The optimal allocation models we consider can be formulated in the standard
optimal control setting [30, 42]. The state variables are $x$ (vegetative
mass) and $y$ (reproductive mass), and the net rate of photosynthesis depends
only on the vegetative mass given by the production function $P(x)$. The
control variable $u$ is the fraction of the net production reinvested into
vegetative growth, and the remainder $1-u$ is invested in reproduction. The
starting time $t=0$ is the time when the plants shift from maternal seed
resources to photosynthesis for growth. We assume that plants’ lifetimes are
identically distributed, but not necessarily independent. Unfortunately, using
something like a Gaussian distribution makes the problem analytically
intractable. As a compromise, we split the lifetime into a “safe” period, with
fixed length $T_{0}$, and the “volatile” period of maximal length $T$, which
is equally likely to end at any time after $T_{0}$. The state equations are as
follows:

$\displaystyle\begin{cases}\dot{x}=uP(x),\ x(0)=x_{0}>0\\\ \dot{y}=(1-u)P(x),\
y(0)=y_{0}\geq 0;\end{cases}$ (1) $\displaystyle 0\leq u\leq 1,\ \ \ \ \
t\in[0,T_{0}+T].$

The payoff functional is the expected value of the fitness function $L(y)$
over possible lifetimes (we drop the $\frac{1}{T}$ factor that does not affect
optimization):

$J[u](x_{0},y_{0},T_{0},T)=\int_{T_{0}}^{T_{0}+T}L(y)dt\to\max\,.$ (2)

We generally assume that $P(x)$ is positive, and both $P(x)$ and $L(y)$ are
twice differentiable and monotone increasing for positive $x$ and $y$. Other
conditions will be explicitly stated when necessary. Under these conditions
the existence of optimal control follows from the Filippov’s theorem [42,
p.132].

To state the results more compactly we use a shorthand notation for different
arcs of optimal allocation schedules. The letters below are V for vegetative,
M for mixed, and R for reproductive growth.

######  Definition 1.

We label an arc of an optimal schedule V when $u=1$, M when $0<u<1$, and R
when $u=0$, respectively. A schedule is said to be of type VMR, etc., when the
corresponding arcs are arranged from $t=0$ to $t=T$ in the indicated order
left to right.

In principle, a schedule can include any sequence of arc types in any order,
but their variety turns out to be very limited for optimal schedules. The
following theorem characterizes them under very general assumptions on
production and fitness functions.

######  Theorem 1.

Let $\ln P$ be concave and $L$ be strictly monotone increasing. Then

(i) If $L$ is convex and $L(y)\xrightarrow[y\to\infty]{}\infty$ then any
optimal schedule has at most one switch, and is bang-bang, i.e. R or VR;

(ii) If $\ln L^{\prime}$ is strictly convex then any optimal schedule has no
more than three switches, and one of the following schedule types: R, VR, VMR,
RMR, VRMR. Moreover, if $T_{0}$ is not too short (e.g.
$T_{0}P^{\prime}(x(T_{0}))\geq 1$) and $y(0)=0$ then only VR and VRMR
schedules are possible, and the first switch occurs during the safe period.

If $T_{0}$ is very short the plant may have to start reproducing from the
start, which is not biologically feasible. The value of $P^{\prime}(x(T_{0}))$
is typically easy to estimate from below by integrating $\dot{x}=P(x)$ on
$[0,T_{0}]$. Since $y(0)=0$ means that the plant has no reproductive mass at
the time of establishment, only VR and VRMR schedules have biological
relevance when $L$ is as in Theorem 1(ii). The first switching time $t_{m}>0$
is the time when the plant starts reproducing for the first time, and can be
identified with the age of maturity. For realistic parameter values it is
often close to $T_{0}$. When the schedule is not bang-bang, i.e. when it is
VRMR, the first reproductive phase is often short and seems to be an artifact
of our sharp separation between safe and volatile periods. With a more
continuous distribution of lifetimes one can expect a direct transition from
vegetative to mixed growth followed by purely reproductive growth.

To interpret the conditions on $L$ more concretely let us apply them to the
power functions $L(y)=\frac{y^{\alpha}}{\alpha}$. They are convex for
$\alpha\geq 1$ corresponding to “gambler’s” behavior, and $\ln
L^{\prime}(y)=(\alpha-1)\ln y$ is strictly convex for $\alpha<1$. In other
words, the possibility of mixed growth only exists for the more conservative
indices $\alpha<1$. In fact, if the volatile period is long enough then mixed
growth always occurs for some initial values. A quantitative threshold can be
estimated in terms of two constants $q_{c},\tau_{c}$ that depend only on
$\alpha$, which we call the cutoff constants. They are easiest to interpret
when the production function is linear with unit rate, $P(x)=x$. Then
$\tau_{c}$ is the time-to-go at the end of the mixed period, i.e. (maximal)
time remaining until the end of plant’s lifetime, and $q_{c}$ is the ratio of
reproductive to vegetative mass at that time. The constants are easy to find
to any precision by solving a system of two non-linear equations, see Section
4.2. For $\alpha=1,\pm\frac{1}{2}$ and $-2$ they can be found analytically,
and both grow to $\infty$ when $\alpha\to-\infty$. Some representative values
are given in Table 1, for $L(y)=\ln y$ they were found in [28].

######  Theorem 2.

Let $P(x)$ be concave and $L(y)=\frac{y^{\alpha}}{\alpha}$ with $\alpha<1$. If
$T>\frac{q_{c}+\tau_{c}}{P^{\prime}(x(T_{0}))}$ then any optimal schedule with
$y(0)=0$ has a mixed growth arc. In particular, when $P(x)=kx^{\beta}$ it
suffices that
$T>\frac{q_{c}+\tau_{c}}{k\beta}\big{(}x_{0}^{1-\beta}+(1-\beta)kT_{0}\big{)}$.

$\alpha$ | $1$ | 1/2 | $0$ | $-1/2$ | $-2$
---|---|---|---|---|---
$q_{c}$ | $0$ | 1/6 | $0.558$ | 1 | $1+\sqrt{2}$
$\tau_{c}$ | $2$ | 5/2 | $2.793$ | 3 | $2+\sqrt{2}$
Table 1: Cutoff constants that determine occurrence of mixed growth.

We should point out that unlike the condition on $T_{0}$, which is always
satisfied in realistic circumstances, $T$ may well not be “long enough”. In
fact, the minimal length explicitly depends on $\alpha$, and mixed phase only
occurs for a limited range of its values, other parameters being equal.

Much more detailed results are obtained for linear production functions
$P(x)=kx$, when the model can be solved (almost) analytically, see Section 5.
The optimal growth trajectories are fully characterized by the evolution of
mass ratio $q=\frac{y}{x}$, and the full picture of their dependence on
initial values and model parameters can be visualized on a 2D diagram. This
generalizes the results of [28] for $L(y)=\ln y$.

###  1.3 Discussion

Since Cohen’s work in 1960s, geometric mean maximization was associated with
risk aversion and conservatism, while arithmetic mean maximization with more
risk taking. It is easy to see why. In a large well-mixed population when
lifetimes are uncorrelated we get a large number of essentially independent
trials for the same allele type. Even if many plants die without offspring,
survivors leave more than enough to compensate. On the other hand, when the
correlation is perfect population size does not matter much because all plants
die almost at once (e.g. due to water shortage or a cold spell). A vivid
illustration of this point is that geometric mean becomes $0$ if even one of
the terms is $0$. But in general, seasonal and individual adult mortality
factors (like grazing by animals or microsite environment) combine and the
correlation is neither perfect nor zero.

We argued that some non-linear effects of density dependence (such as lifetime
variability and seedling mortality) can be represented by introducing a non-
linear fitness function $L(y)$, whose mathematical expectation is maximized.
Its role is analogous to the role of utility function in economics, but only
linear and logarithmic functions were previously considered in the ecological
literature. When reproductive mass is proportional to the effective seed yield
$L(y)=\frac{y^{\alpha}}{\alpha}$ was shown to be a natural choice, and the
fitness index $\alpha$ being directly related to the correlation coefficient
$\rho$ between lifetime distributions, $\alpha=1-\rho$. By varying $\alpha$
from $0$ to $1$ (correspondingly, $\rho$ from $1$ to $0$) our model traces
changes in the optimal schedule. If juvenile mortality saturates effective
seed yields, negative values of $\alpha$ become meaningful as well, while seed
predation and low disperser efficiency at low seed yields may lead to
“gambler’s” fitness with $\alpha>1$.

Two observations stand out. First, the time of maturity, when reproduction
starts, grows with $\alpha$, moving from within the safe period for small
values to within the volatile period for $\alpha\geq 1$. Second, while on both
ends optimal schedules are bang-bang, the transition is mediated by schedules
with mixed growth for a wide range of intermediate values. On the upper end
the mixed period disappears at $\alpha=1$ for any values of other parameters,
and does not appear for “gambler’s” fitness $\alpha>1$, but on the lower end
the threshold value depends on the initial plant size and the lifetime
distribution.

(a) (b) (c)
(d) (e) (f)

Figure 2:  Optimal schedules (graphs of $u(t)$) for $P(x)=0.1x^{0.8}$ with
$T_{0}=40$ and $T=80$ days, $x_{0}=2.5$ mg, $y_{0}=0$, and (a) $\alpha=-4$ (b)
$\alpha=-3$ (c) $\alpha=-0.5$ (d) $\alpha=0.1$ (e) $\alpha=0.7$ (f)
$\alpha=1.5$.

Figure 2 shows a progression of optimal schedules as $\alpha$ increases for a
choice of parameter values. While thresholds and durations vary with model
parameters the pattern is typical. The duration of mixed period increases with
$\alpha$, and its nature shifts from reproductive growth dominating at smaller
values to vegetative growth dominating as $\alpha$ approaches $1$. For the
chosen parameter values it first appears at $\alpha\approx-3.5$, and at
$\alpha=1$ it merges with the initial period of vegetative growth. No mixed
growth occurs for $\alpha>1$, in accordance with Theorem 1(i). In the course
of mixed growth period the fraction of net production allocated to
reproduction increases with time. The intermediate reproductive period, that
precedes mixed growth, is barely noticeable for a range of $\alpha$, in this
case $0<\alpha<1$, and is likely an artifact of the sharp separation between
safe and volatile periods in our model. In this range the switch to mixed
growth is very close to the beginning of the safe period.

For a continuous distribution of lifetimes one can still expect that the time
of maturity moves from the lower tail of the lifetime distribution towards the
average as $\alpha$ increases, and the initial vegetative growth transitions
directly into mixed growth. It then follows a balanced path where the marginal
values of vegetative and reproductive output are equal. The model also
predicts that mixed growth is often followed by an opportunistic reproductive
spurt at the end of plant’s lifetime.

Mechanisms that support mixed growth in deterministic models include
increasing the number of allocation sectors [25, 34], introducing size-
dependent mortality and/or fecundity [26, 36], and size-dependent constraints
on the size of reproductive organs [24, 49]. As we saw, fluctuating
environments lead to mixed growth being optimal under certain conditions even
without any additional mechanisms. In continuous time models this phenomenon
was first discovered by King and Roughgarden [28] for geometric mean
maximizers and linear production. Our results extend the range where it is to
be expected, and clarify limiting conditions on its occurrence. It is robust
over a range of production functions with diminishing returns. It is ruled out
for linear and convex fitness functions, as well as saturating ones that
diminish selective advantage returns on reproduction too quickly. This last
threshold is sensitive to the initial mass and lifetime distribution of plants
in addition to the fitness index.

The structure of optimal trajectories $x(t),y(t)$ has a nice geometric
description in terms of singular and switching surfaces that collect balanced
growth paths and switching points between vegetative and reproductive growth
periods, respectively. For linear production functions optimal dynamics
depends only on the ratio of vegetative and reproductive masses, and can be
presented on a two-dimensional diagram (Figures 9 and 12), that visually
manifests many features of the model predictions valid also in the general
case. For $L(y)=\ln y$ such a diagram was presented already in [28]. Non-
linear production shifts the onset of various stages depending on the initial
values, but the overall pattern remains the same. As can be expected, mixed
growth period lasts longer when the production function is made proportionally
bigger, or less concave, because both provide greater opportunities for
leveraging larger plant size when generating seeds.

In summary, our results suggest that when the length of lifetime fluctuates
bang-bang behavior predicted by deterministic models cannot be taken as the
baseline against which other mechanisms that induce mixed growth are assessed.
The interaction of risk-averse bet-hedging with such mechanisms, as well as
with seed bank effects, is of great interest for further study.

The rest of the paper is mostly mathematical and is organized as follows. In
Sections 2-5 we study the model for the volatile period only, which is most
challenging mathematically. In particular, in Section 2 we introduce the
marginal values of vegetative and reproductive masses, and characterize
optimal schedule types, as needed to prove Theorem 1. In Section 3 we
introduce geometric objects that determine which schedule type is selected for
given initial values, and establish some of their properties. In Section 4 we
specialize to power fitness functions, and prove the main ingredients for
Theorem 2. Section 5 presents a detailed analytic solution for the case of
linear production. Finally, in Section 6 we incorporate the volatile period
results into the full lifetime model.

## List of symbols

$x,y$ | vegetative, reproductive mass, 1.2 |  | $\lambda,\mu$ | costate variables for $x$, $y$, 3.1
---|---|---|---|---
$T_{0},T$ | safe, volatile period, 1.2 |  | $\xi$ | adjusted vegetative mass, 3.1
$L(y)$ | fitness function, 1.1 |  | $\tau$ | adjusted time-to-go, 3.1
$P(x)$ | production function, 1.2 |  |  $q$, $r$ | adjusted mass ratios, 4
$\alpha$ | fitness index, 1.1 |  |  $q_{c}$, $r_{c}$, $\tau_{c}$ | cutoff constants, 4; 4.2
$\textrm{Sw}(\tau,\xi,y)$ | switching function, 3.1 |  | $t_{s}$ | VR switching time, 1.2; 6
$\textrm{Cut}(\tau,\xi,y)$ | cutoff function, 3.1 |  | $t_{m}$ | time of maturity, 3.1

##  2 Optimal schedule types for the volatile period

Most of the paper studies growth during the volatile period, which is more
challenging mathematically, and the full lifetime model is considered in
Section 6. Until then, we take the volatile period to start at $t=0$, but it
is important to keep in mind that the plant is likely to be mature already at
the starting point.

In the following we use the standard notation and terminology of the optimal
control theory [30, 42]. The payoff functional for the volatile period
simplifies to:

$J[u](x_{0},y_{0},T)=\int_{0}^{T}L(y)dt\to\max\,.$ (3)

The costate (adjoint) variables are denoted $\lambda$ and $\mu$, and represent
the marginal values of increasing $x$ and $y$, respectively. In other words,
$\lambda(t)$ is the increase in the payoff from a unit increase in the
vegetative mass at time $t$, and similarly for $\mu(t)$ [30]. The Hamiltonian
of the model is

$H(x,y,\lambda,\mu,u)=L(y)+\big{(}\mu+(\lambda-\mu)u\big{)}P(x).$ (4)

The costate equations and the transversality conditions at the terminal time
$T$ are given by Pontryagin’s theorem:

$\begin{cases}\dot{\lambda}=-\big{(}\mu+(\lambda-\mu)u\big{)}P^{\prime}(x),\ \
\lambda(T)=0\\\ \dot{\mu}=-L^{\prime}(y),\ \ \mu(T)=0.\\\ \end{cases}$ (5)

Since our system is autonomous, it satisfies the stationarity condition, i.e.
the Hamiltonian remains constant along any extremal trajectory. The
transversality conditions and (4) imply that this constant is equal to
$L(y(T))$. Since the state equations and the payoff functional are linear in
control we are dealing with a control-affine system, a class which is well
known to feature both bang-bang and singular controls [30].

According to the Pontryagin maximum principle, the optimal control is obtained
by maximizing the Hamiltonian as a function of $u$ at each time $t$. Since
$\frac{\partial H}{\partial u}=(\lambda-\mu)P(x)$ does not depend on the
control explicitly, the maximum of $H$ as a function of $u$ is attained at one
of the endpoints, $0$ or $1$, as long as $\frac{\partial H}{\partial u}\neq
0$. Since $P(x)>0$ for positive $x$ we have:

$u(t)=\begin{cases}1,\ \lambda(t)>\mu(t)\\\ 0,\
\lambda(t)<\mu(t)\,.\end{cases}$ (6)

Therefore, singular controls $0<u<1$ corresponding to mixed growth are only
possible when $\lambda(t)=\mu(t)$ over an entire interval, in which case
$\dot{\lambda}(t)=\dot{\mu}(t)$. It follows from the costate equations that
$\lambda(t)$ and $\mu(t)$ are continuously differentiable as functions of $t$,
as long as the control is piecewise continuous. Moreover,
$\mu=\frac{L^{\prime}(y)}{P^{\prime}(x)}$ on any singular arc, and
differentiating this identity while taking into account the state equations
yields a feedback expression for the singular control:

$u=\frac{L^{\prime\prime}(y)+L^{\prime}(y)\frac{P^{\prime}(x)}{P(x)}}{L^{\prime\prime}(y)+L^{\prime}(y)\frac{P^{\prime\prime}(x)}{P^{\prime}(x)}}.$
(7)

While equations (1) with this expression for singular $u$ can not usually be
integrated explicitly, some useful information can be extracted from the
stationarity condition. On a singular trajectory the Hamiltonian reduces to

$H(x,y)=L(y)+\frac{P(x)}{P^{\prime}(x)}L^{\prime}(y),$ (8)

and remains constant, which gives an algebraic relation between $x$ and $y$
that holds during periods of mixed growth.

We will now determine possible switching structures of optimal control in the
allocation problem (1)-(3). Although a priori any sequence of arc types in any
order is possible, it turns out that the possibilities in our case are quite
limited. This is largely because the costate equations guarantee that both
$\lambda(t)$ and $\mu(t)$ are strictly monotone decreasing. We summarize the
result in the following theorem. Note that $\ln P$ is concave when $P$ is
since $P$ and $P^{\prime}$ are positive.

######  Theorem 3.

Let $P(x)$ and $L(y)$ be twice differentiable and $P(x)$ be positive for
$x,y>0$. Then we have for the optimal trajectories of problem (1), (3):

(i) If $L(y)$ is strictly monotone increasing, i.e. $L^{\prime}(y)>0$, then
the final arc of an optimal trajectory is an R arc;

(ii) If $\ln P$ is concave, and an optimal trajectory has a V arc, this arc is
initial;

(iii) If $\ln L^{\prime}$ is strictly convex, and an optimal trajectory has an
R arc preceding an M arc, this R arc is initial;

(iv) If $P(x)$ and $L(y)$ meet all the conditions of (i)-(iii) then any
optimal trajectory has no more than two switches, and one of the following
control structures: R, MR, VR, RMR, VMR.

###### Proof.

(i) From the transversality conditions, $\mu(T)=\lambda(T)=0$, and the adjoint
equations (5), we have that $\dot{\lambda}(T)=0$, and
$\dot{\mu}=-L^{\prime}(y(T))<0$. Therefore, in a small left neighborhood of
$T$ we have $\lambda(t)<\mu(t)$, and $u=0$ by the maximum principle.

(ii) By contradiction, supposed that there is a switch to the V arc at time
$t_{s}>0$. We have that $\lambda(t_{s})$ = $\mu(t_{s})$ for the switching time
$t_{s}$. As the arc before $t_{s}$ is V, we know that $\lambda$ must have been
greater than $\mu$ for $t$ close to $t_{s}$, so that
$\dot{\lambda}(t_{s})\geq\dot{\mu}(t_{s})$. Taking the derivative of
$\eqref{Adj}$ when $u=1$, we get
$\ddot{\lambda}=(P^{\prime}(x)^{2}-P^{\prime\prime}(x)P(x))\lambda$. Since
$(\ln
P)^{\prime\prime}=\frac{-P^{\prime}(x)^{2}+P^{\prime\prime}(x)P(x)}{P^{2}(x)}<0$
we have $\lambda>0$, and
$\ddot{\mu}=\frac{d}{dt}\big{(}-L^{\prime}(y)\big{)}=0$, since $y$ is
constant. Thus, if $\lambda>\mu$, $\mu$ and $\lambda$ cannot intersect again
on $[0,t_{s})$, as claimed.

Figure 3:  Marginal values on $[0,t]$ when $\mu(t)>\lambda(t)$ for $t<t_{s}$.
The values are equal for $t\geq t_{s}$ (purple curve), so $t_{s}$ is an RM
switching time.

(iii) Let $t_{s}0$ be the time of the RM switch. We have that
$\lambda(t_{s})=\mu(t_{s})=\dot{\lambda}(t_{s})=\dot{\mu}(t_{s})=-L^{\prime}(y_{s})$.
As the arc before $t_{s}$ is R, we know that $\lambda<\mu$ on a small interval
immediately before $t_{s}$. Furthermore, since $\lambda(t_{s})=\mu(t_{s})$ and
$\dot{\lambda}(t_{s})=\dot{\mu}(t_{s})$, we know
$\ddot{\lambda}(t_{s})\leq\ddot{\mu}(t_{s})$. Solving $\eqref{Adj}$ with
$u=0$, we reduce this inequality to:

$\frac{P^{\prime}(x(t_{s}))}{P(x(t_{s}))}\leq-\frac{L^{\prime\prime}(y(t_{s}))}{L^{\prime}(y(t_{s}))}<-\frac{L^{\prime\prime}(y(t))}{L^{\prime}(y(t))}.$

The last inequality follows from $y(t)$ strictly increasing and
$\frac{L^{\prime\prime}}{L^{\prime}}$ strictly decreasing for $t<t_{s}$ since
$\frac{L^{\prime\prime}}{L^{\prime}}=(\ln(L^{\prime}))^{\prime}$. Therefore,
$\ddot{\lambda}(t)\leq\ddot{\mu}(t)$, and hence $\lambda(t)<\mu(t)$ for
$t<t_{s}$, see Figure 3. Thus, there are no earlier switches, and the R arc is
initial.

(iv) By (i), the last arc is always R. If there is no switching we have a pure
R trajectory, and if it is preceded immediately by a V arc then the latter is
initial by (ii), so we have a VR trajectory. If it is immediately preceded by
an M arc, then either M is initial, and we have an MR trajectory, or the M arc
is preceded by a V or R arc. Either way, those are initial by Lemmas (ii),
(iii), respectively. ∎

Already for $P(x)=x$ and $L(y)=y^{\alpha}$ with any $\alpha<1$ there are
values of $x_{0},y_{0}$ and $T$ that instantiate each of the five structures
(Section 5). On the other hand, for $\alpha\geq 1$ only bang-bang, i.e. R and
VR, structures occur (Corollary 1). The next question is how to tell which
structure is optimal for a given set of parameters.

##  3 Geometry of optimal trajectories

It turns out that the switching structure of an optimal schedule is determined
by the location of the initial point relative to a geometric configuration,
two surfaces and the curve of their intersection, see Figure 4. We describe
the elements of this configuration in this section, and sketch how it
determines optimal strategies at the end of 3.3. For linear production
functions the configuration collapses to a two-dimensional one described in
Section 5.

###  3.1 Switching surface

Let us start by looking at the VR trajectories. The marginal values are both
equal (to zero) at the terminal time $T$. For the final R arc to join to the V
arc at an earlier time, they must become equal again at that time. This
motivates the following definition.

######  Definition 2.

We say that a point $(x,y,t)$ is on the switching surface if $t<T$, and the
costate variables computed backwards in time along a pure R extremal
trajectory starting at it, are equal.

The trajectory in the definition is uniquely determined. Denote $t_{s}$ the
switching time, and $x_{s},y_{s}$ the values at the switch. Since $u=0$ after
the switch we have that $x(t)=x_{s}$ is constant, and
$y(t)=y_{s}+P(x_{s})(t-t_{s})$. It is convenient to use $y$ as the integration
variable when integrating the costate equations (5) backwards from $T$ to
$t_{s}$. Denoting $y_{T}=y(T)$, we obtain:

$\mu(y)=\frac{L(y_{T})-L(y)}{P(x_{s})},\quad\lambda(y)=\frac{P^{\prime}(x_{s})}{P(x_{s})^{2}}\int^{y_{T}}_{y}L(y_{T})-L(y)\,dy.$
(9)

The difference $T-t_{s}$ is often called time-to-go in control theory. Since
$\lambda(t)=\mu(t)$ at the VR junction, we obtain a relation between the
switching values $t_{s},x_{s}$ and $y_{s}$ satisfied at such a junction. To
express it, it is convenient to introduce new variables:

$\displaystyle\xi:=P(x)/P^{\prime}(x),\ \ \ \tau:=P^{\prime}(x)(T-t),\ \ \
y_{T}=y+\xi\tau,$ (10)

and define the switching function:

$\displaystyle\textrm{Sw}(\tau,\xi,y)$
$\displaystyle:=L(y_{T})(y_{T}-y)-\int_{y}^{y_{T}}L(\zeta)d\zeta-\xi\big{(}L(y_{T})-L(y)\big{)}$
(11) $\displaystyle=\int_{y}^{y_{T}}\zeta
L^{\prime}(\zeta)d\zeta-(\xi+y)\big{(}L(y_{T})-L(y)\big{)}\,.$

Note that in $\tau$, $\xi$, $y$ coordinates the equation of the switching
surface does not depend on $P(x)$. One can see that

$\textrm{Sw}(\tau_{s},\xi_{s},y_{s})=\frac{1}{\xi_{s}P(x_{s})}\big{(}\lambda(y_{s})-\mu(y_{s})\big{)},$
(12)

i.e. the switching surface is given by the equation $\textrm{Sw}=0$ with
$\tau>0$. This equation allows one to find VR trajectories (almost)
analytically. On the V arc we have $\dot{x}=P(x)$ and $\dot{y}=0$. This gives
$y_{s}=y_{0}$ and a system of two equations to determine $t_{s}$ and $x_{s}$:

$t_{s}=\int\limits_{x_{0}}^{x_{s}}\frac{dx}{P(x)};\ \ \
\textrm{Sw}\left(P^{\prime}(x_{s})(T-t_{s}),\frac{P(x_{s})}{P^{\prime}(x_{s})},y_{s}\right)=0.$
(13)

Once these values are found, the VR trajectory can be easily integrated, and
satisfies the necessary conditions of the Pontryagin’s theorem. Of course,
these trajectories are not always optimal.

Let us prove some simple properties of the switching surface.

######  Lemma 1.

A trajectory that passes through a point above the switching surface (i.e.
$\textrm{Sw}>0$), and is R afterwards, is suboptimal.

###### Proof.

Without loss of generality, we may assume that at our point $t=0$. Let
$J(x,y,T)$ denote the payoff computed along an R trajectory starting from
$(x,y)$ over $[0,T]$. To show that such a trajectory is suboptimal we will
produce an admissible trajectory with higher payoff. Consider a trajectory
with $u=1$ on $[0,\epsilon]$, and $u=0$ on $[\epsilon,T]$. Denoting
$x_{\epsilon}=x(t_{\epsilon})$ we can express the new payoff as:

$J_{\epsilon}=L(y)\epsilon+J(x(\epsilon),y,T-\epsilon).$

Now $\frac{dx(\epsilon)}{d\epsilon}\Big{|}_{\epsilon=0}=P(x)$, $\frac{\partial
J}{\partial T}=L(y_{T})$, and $\frac{\partial J}{\partial x}=\lambda(0)$
computed along the R trajectory (by the standard properties of marginal values
[30, 2.2]). Taking into account (9), we compute

$\frac{dJ_{\epsilon}}{d\epsilon}\Big{|}_{\epsilon=0}=P(x)\big{(}\lambda(0)-\mu(0)\big{)}.$

Since $(x,y,0)$ is above the switching surface, (12) implies that
$\lambda(0)>\mu(0)$, and the derivative is positive. But then
$J_{\epsilon}>J_{0}$ for small $\epsilon>0$, as required. ∎

Next, we will show that for monotone non-increasing functions $L(y)$ the
switching surface does not come too close to the plane $t=T$, i.e. the time-
to-go at the VR junction cannot be too small. This is plausible intuitively
since the plant needs enough time to convert accumulated $x$ into $y$. The
inequality we will prove is sharp since, for $L(y)=y$ and $P(x)=x$, the
switching surface is the plane $T-t=2$ (VR switches in this case happen with a
constant time-to-go as follows from (21) with $\alpha=1$, for example).

######  Theorem 4.

Suppose that $L^{\prime}(y)$ is monotone non-increasing. Then the time-to-go
on the switching surface (at the junction with the final R arc) satisfies
$T-t_{s}\geq\frac{2}{P^{\prime}(x_{s})}$. This inequality is strict if
$L^{\prime}(y)$ is strictly decreasing.

###### Proof.

It will be convenient to use the ratio variable $r:=\frac{\xi}{y}$. We will
first transform the second equation in (11). Writing
$L(y_{T})-L(y)=\int^{y_{T}}_{y}L^{\prime}(\zeta)d\zeta$, and changing to a new
variable $s=\frac{\zeta}{y}$ in both integrals, we obtain:

$\textrm{Sw}(\tau,\xi,y)=y^{2}L^{\prime}(y)\left(\int_{1}^{1+r\tau}\\!\\!\\!\\!\\!\\!sL^{\prime}(ys)ds-(1+r)\int_{1}^{1+r\tau}\\!\\!\\!\\!\\!\\!L^{\prime}(ys)ds\right)\,.$
(14)

Next, we will need the following elementary inequality for positive monotone
non-increasing functions $l(s)$:

$\int_{1}^{R}\\!\\!\\!\\!s\,l(s)ds\leq\frac{R+1}{2}\int_{1}^{R}\\!\\!\\!\\!l(s)ds$

Indeed, this is equivalent to
$\displaystyle{\int_{1}^{R}\\!\\!\\!\left(s-\frac{R+1}{2}\right)l(s)ds\leq
0}$, which holds as the difference $s-\frac{R+1}{2}$ is symmetric with respect
to the middle of the interval $[1,R]$, and larger values of $l(s)$ are on the
side where it is negative. The equality is only attained when
$l=\textrm{const}$.

Therefore, taking $l(s)=L^{\prime}(ys)$ in (14) we have that for the points on
the switching surface:

$\int_{1}^{1+r\tau}\\!\\!\\!\\!\\!\\!sL^{\prime}(ys)ds=(1+r)\int_{1}^{1+r\tau}\\!\\!\\!\\!\\!\\!L^{\prime}(ys)ds\leq\frac{2+r\tau}{2}\int_{1}^{1+r\tau}\\!\\!\\!\\!\\!\\!L^{\prime}(ys)ds.$

Canceling the integrals and $r$, we obtain $2\leq\tau$. The conclusion now
follows from (10). ∎

###  3.2 Cutoff curve

Not every point on the switching surface can be the VR junction of an optimal
trajectory. When it is, it follows from (6) that $\lambda>\mu$ reverses to
$\lambda<\mu$ there. Therefore, at the junction points not only does
$\lambda=\mu$, but also $\dot{\lambda}\leq\dot{\mu}$.

######  Definition 3.

A point $(x,y,t)$ of the switching surface is said to be on the cutoff curve
when the time derivatives of the costate variables computed along a pure R
extremal trajectory starting at it are equal, $\dot{\lambda}=\dot{\mu}$. If
$\dot{\lambda}<\dot{\mu}$ we say that the point is above the cutoff curve.

Intuitively, the cutoff curve separates the “active” part of the surface,
where the optimal junctions can occur, from the “inactive” part. Optimal VR
junctions can only occur on or above the cutoff curve, hence the name. As we
will see below, this curve serves as a “cutoff” also for singular arcs. The
cutoff condition can be expressed explicitly by using (9), and recalling that
for $\lambda=\mu$ we have:

$\dot{\lambda}=-\mu P^{\prime}(x_{s})=-L^{\prime}(y)=\dot{\mu}.$

Recall from (10) that $y_{T}=y+\xi\tau$, and define the cutoff function:

$\textrm{Cut}(\tau,\xi,y)=L(y_{T})-L(y)-\xi L^{\prime}(y).$ (15)

Note that by the costate equations

$\textrm{Cut}(\tau_{s},\xi_{s},y_{s})=-\frac{1}{\xi_{s}}\big{(}\dot{\lambda}(y_{s})-\dot{\mu}(y_{s})\big{)}.$
(16)

Thus, points on the cutoff curve can be characterized as satisfying
$\textrm{Cut}=0$, and those above it as satisfying $\textrm{Cut}>0$.

As already mentioned, the cutoff curve is also a cutoff in another sense. As
the only place where $\lambda=\mu$ and $\dot{\lambda}=\dot{\mu}$ on the final
R arc, it is also the only place where an optimal junction with singular M
arcs can occur. At these junctions, the singular control (7) must be
terminated and switched to $u=0$. Conversely, if the cutoff curve is empty, or
the singular control on it is inadmissible, optimal trajectories can only be
bang-bang. The next theorem gives sufficient conditions for the cutoff curve
to be either empty or “large”.

######  Theorem 5.

Suppose that $L(y)$ is strictly increasing and
$L(y)\xrightarrow[y\to\infty]{}\infty$. If $L(y)$ is convex then the cutoff
curve is empty. If $\frac{L(y)}{y}\xrightarrow[y\to\infty]{}0$ then for any
$y>0$ there is a point on the cutoff curve with this value of $y$. If,
moreover, $L(y)$ is strictly concave, this point is unique.

###### Proof.

Since $\textrm{Sw}=\textrm{Cut}=0$ on the cutoff curve, we get from (11) and
(15):

$\int_{y}^{y_{T}}\zeta L^{\prime}(\zeta)d\zeta=(\xi+y)\xi L^{\prime}(y)$ (17)

By our assumptions, $L(y)$ has an inverse function $L^{-1}(z)$ defined for
$z>0$, and $L^{-1}(z)\to\infty$ when $z\to\infty$. By making a change of
variables in the integral, $\zeta=L(y)$, taking into account that
$L(y_{T})=L(y)+\xi L^{\prime}(y)$, and denoting $\eta:=\xi L^{\prime}(y)$, we
get the following equation for $\eta$:

$\int^{L(y)+\eta}_{L(y)}\\!\\!L^{-1}(z)dz=\Big{(}\frac{\eta}{L^{\prime}(y)}+y\Big{)}\eta.$
(18)

Let $F(\eta)$ and $G(\eta)$ denote the left and the right hand sides of (18),
respectively. Then we have $F(0)=G(0)=0$, and

$\displaystyle
F^{\prime}(\eta)\Big{|}_{\eta=0}=L^{-1}\big{(}L(y)+\eta\big{)}\Big{|}_{\eta=0}=y;\
\ \ \ \ \ \ \ \ \ \ \ \ \
G^{\prime}(\eta)\Big{|}_{\eta=0}=y+2\frac{\eta}{L^{\prime}(y)}\Big{|}_{\eta=0}=y$
$\displaystyle
F^{\prime\prime}(\eta)\Big{|}_{\eta=0}=(L^{-1})^{\prime}\big{(}L(y)+\eta\big{)}\Big{|}_{\eta=0}=\frac{1}{L^{\prime}(y)};\
\ \ \ \ G^{\prime\prime}(\eta)\Big{|}_{\eta=0}=\frac{2}{L^{\prime}(y)}.$ (19)

Therefore, $F^{\prime\prime}(\eta)<G^{\prime\prime}(\eta)$,
$F^{\prime}(\eta)<G^{\prime}(\eta)$, and $F(\eta)<G(\eta)$ for small $\eta>0$.
When $L(y)$ is convex, $L^{-1}(z)$ is concave, and $(L^{-1})^{\prime}(z)$ is
non-increasing, so $(F-G)^{\prime\prime}(\eta)$ is as well. But then, it
remains negative for all $\eta>0$, and so do $(F-G)^{\prime}$ and $F-G$. Thus,
there are no positive solutions to $F(\eta)=G(\eta)$.

When $\frac{L(y)}{y}\xrightarrow[y\to\infty]{}0$, we have
$\frac{L^{-1}(z)}{z}\xrightarrow[y\to\infty]{}\infty$, and
$\frac{F^{\prime}(\eta)-G^{\prime}(\eta)}{\eta}\xrightarrow[\eta\to\infty]{}\infty$.
Therefore, $F(\eta)-G(\eta)\xrightarrow[\eta\to\infty]{}\infty$, and
$F(\eta)=G(\eta)$ for some $\eta>0$. Setting

$\xi:=\frac{\eta}{L^{\prime}(y)},\ \
\tau:=\frac{L^{\prime}(y)}{\eta}\left(L^{-1}\big{(}L(y)+\eta\big{)}-y\right)$

gives us the desired point on the cutoff curve.

If $L(y)$ is strictly concave, then $L^{-1}(z)$ is strictly convex, and
$(L^{-1})^{\prime}(z)$ is strictly increasing. Then, so is
$(F-G)^{\prime\prime}$. If it stayed negative for $\eta>0$ then so would
$(F-G)^{\prime}$ and $F-G$, contradicting the existence of a solution.
Therefore, $F^{\prime\prime}(\eta_{**})=G^{\prime\prime}(\eta_{**})$ for some
$\eta_{**}>0$. Similarly, $F^{\prime}(\eta_{*})=G^{\prime}(\eta_{*})$ for some
$\eta_{*}\geq\eta_{**}$, and the difference is positive for $\eta>\eta_{*}$.
Thus, $F-G$ is strictly increasing for $\eta>\eta_{*}$, and since it is
negative for $\eta<\eta_{*}$, it must be zero at a single positive $\eta$. ∎

######  Corollary 1.

If $L(y)$ is strictly increasing, convex, and
$L(y)\xrightarrow[y\to\infty]{}\infty$, then the optimal control is always
bang-bang, and the optimal trajectories are either R or VR.

This situation occurs in the case $L(y)=y$, and for $L(y)=y^{\alpha}$ with
$\alpha\geq 1$ generally. We see, however, that it is borderline to the range
of gambler’s fitness functions, while the more common concave fitness
functions have non-empty cutoff curves. The increase of $L(y)$ to $\infty$ is
a technical condition, one can modify the proof to cover functions with a
horizontal asymptote, like $L(y)=\frac{y^{\alpha}}{\alpha}$ for $\alpha<0$,
see Lemma 3.

###  3.3 Singular surface

We are now in a position to define the final piece of our geometric
configuration, the singular surface. When we substitute the singular control
law (7) into the state equations (1) we obtain integral curves passing through
every point of state space. However, these curves can only produce singular
arcs of optimal trajectories if they pass through a point on the cutoff curve,
where these arcs can be joined to the final R arc. This motivates the next
definition.

######  Definition 4.

We call the singular surface the set foliated by the integral curves of the
system (1) with the singular control (7) (singular curves, for short) passing
through points on the cutoff curve.

By definition, the singular surface is empty when the cutoff curve is empty,
and, when it is non-empty, it intersects the switching surface along the
cutoff curve, see Figure 4. The part of the singular surface past the cutoff
curve is “inactive”, i.e. singular arcs of optimal trajectories can not pass
through it there.

While the non-emptiness of the cutoff curve is a necessary condition for the
existence of singular arcs, it is not sufficient. A priori, the singular
curves may not even be locally optimal, or the singular controls on them may
not be admissible, that is, remain within the bounds $0\leq u\leq 1$.

Figure 4:  Schematic depiction of the switching surface (blue), the cutoff
curve (black), and the singular surface (red) in $\xi$-$q$-$\tau$ coordinates,
$\xi=\frac{P(x)}{P^{\prime}(x)}$ and $q=\frac{y}{\xi}$. The lined parts are
inactive.

We will now sketch how the geometric configuration, as depicted on Figure 4,
determines the nature of optimal schedules. To interpret this biologically,
recall that our model only covers the volatile part of the season, with
uncertain duration, so that the plant is likely to be mature already at the
starting point. Note that the direction of $\tau$ is opposite to the direction
of time $t$, so the progression of time corresponds to moving to the left in
the figure.

Trajectories that start under the singular surface, or to the left of the
switching surface, move up along an R arc until they meet the active part of
the singular surface, or reach $\tau=0$ ($t=T$). In the latter case, they are
pure R, in the former, they switch to a singular M arc, and follow it until
meeting the cutoff curve. When there, they switch to a final R arc, making
them RMR. Trajectories that start above both surfaces move down along a V arc
until meeting one of them (or both, on the cutoff curve). At the switching
surface they switch to a final R arc (VR), and at the singular surface they
first switch to a singular M arc, and then to a final R arc at the cutoff
curve (VMR). MR trajectories must start on the active part of the singular
surface. Of course, this intuitive picture depends on the configuration being
roughly as depicted. In the next section, we will confirm that this is indeed
the case for the power fitness functions, and add some analytic
characterizations.

##  4 Power fitness functions

In this section we apply geometric considerations of the previous section to
the optimal allocation problem (1)-(3) with the power fitness functions
$L(y)=\frac{y^{\alpha}}{\alpha}$. The $\frac{1}{\alpha}$ is added for
technical convenience and does not affect optimization. It simplifies some
formulas, provides the correct sign for $\alpha<0$, and ensures continuity
with $L(y)=\ln y$ for $\alpha\to 0$, as
$\frac{y^{\alpha}-1}{\alpha}\xrightarrow[\alpha\to 0]{}\ln y$.

The geometric configuration we introduced simplifies considerably when $L(y)$
is a power function or logarithm. The underlying reason is that those are
exactly the functions for which the fitness density $L^{\prime}(y)$ is
multiplicative, i.e. $L^{\prime}(ys)=L^{\prime}(y)L^{\prime}(s)$, and the
scaling factor $l_{y}(s):=\frac{L^{\prime}(ys)}{L^{\prime}(y)}$ does not
depend on $y$. As a result, a prominent role is played by the ratio
$r:=\frac{\xi}{y}$ and its reciprocal $q:=\frac{y}{\xi}$. For power production
functions $P(x)=kx^{\beta}$ we have $\xi=\frac{x}{\beta}$, so these ratios are
proportional to the ratios of vegetative and reproductive masses. We will
refer to $q$ as the mass ratio in general. Much of the geometry can be
described in two-dimensional terms, using $q$,$\tau$ or $r$,$\tau$ as
coordinates.

###  4.1 Switching curve

Let us rewrite the switching equation $\textrm{Sw}=0$ in terms of $l_{y}(s)$:

$\int_{1}^{1+r\tau}\\!\\!\\!\\!\\!\\!sl_{y}(s)ds=(1+r)\int_{1}^{1+r\tau}\\!\\!\\!\\!\\!\\!l_{y}(s)ds;$
(20)

Since $l_{y}(s)=l(s)$, the switching surface is a cylinder over a curve in the
$r$-$\tau$ plane given by (20). For $\alpha\neq 0,-1$ it simplifies to111For
$\alpha=0$ this is replaced by $\tau=\left(1+\frac{1}{r}\right)\ln(1+r\tau)$,
and for $\alpha=-1$ by $\tau=\frac{1}{r(1+r)}(1+r\tau)\ln(1+r\tau)$.

$\tau=\left(1+\alpha+\frac{1}{r}\right)\frac{1-(1+r\tau)^{-\alpha}}{\alpha}.$
(21)

We will call the curve given by this equation the switching curve. For
$\alpha=1$ it reduces to $\tau=2$, i.e. the switching curve is a vertical
line. To determine its shape in general, it is convenient to rewrite (21) in a
parametric form using $z=r\tau$ as the parameter (we leave the $\alpha=0,-1$
cases to the reader):

$\displaystyle r(z)$
$\displaystyle=\frac{1}{1+\alpha}\frac{(1+z)^{-\alpha}-1-\alpha
z}{(1+z)^{-\alpha}-1};$ (22) $\displaystyle\tau(z)$
$\displaystyle=(1+\alpha)\,\frac{z\left((1+z)^{-\alpha}-1\right)}{(1+z)^{-\alpha}-1-\alpha
z}.$ (23)

One can prove the following by direct calculations.

######  Lemma 2.

The switching curve is the graph of a positive, monotone decreasing function
$q(\tau)=\frac{1}{r(\tau)}$ on $(2,\infty)$. It has a vertical asymptote at
$\tau=2$, and a horizontal asymptote when $\tau\to\infty$ (and $z\to\infty$).
The latter is $q=0$ for $-1<\alpha<1$ (and the special cases $\alpha=0,-1$),
and $q=-\alpha-1$ for $\alpha<-1$.

For $\alpha=-2$ one can derive an explicit equation,
$q(\tau)=\frac{\tau}{\tau-2}$, see Figure 8.

###  4.2 Cutoff point

In terms of $l_{y}(s)$, the cutoff equation $\textrm{Cut}=0$ becomes

$r=\int_{1}^{1+r\tau}\\!\\!\\!\\!\\!\\!l_{y}(s)ds.$ (24)

Accordingly, the cutoff curve is (generically) a collection of straight lines
over points solving the system (20),(24). We will show when there is, in fact,
such a point which we will call the cutoff point.

######  Lemma 3.

For $\alpha\geq 1$ there are no cutoff points. For $\alpha<1$ there is a
unique point $(r_{c},\tau_{c})$, with $\tau_{c}\geq 2$,
$r_{c}\geq\frac{1}{1-\alpha}$. For $\alpha<0$, the unique point also satisfies
$r_{c}<-\frac{1}{\alpha}$.

###### Proof.

For $\alpha\neq 0,-1$ the system (20),(24) simplifies to222For $\alpha=0$ the
second equation is replaced by $r=\ln(1+r\tau)$, and for $\alpha=-1$ by
$r(1+r)=\ln\tau$.

$\displaystyle\tau$ $\displaystyle=1+\frac{r}{1+\alpha r};$ $\displaystyle
1+\alpha r$ $\displaystyle=(1+r\tau)^{\alpha}.$ (25)

It follows from the second equation that $1+\alpha r\geq 0$ for any positive
solution, and the system reduces to a single equation for $r$:

$f_{\alpha}(r):=(1+\alpha
r)^{\frac{1}{\alpha}}-\left(1+r+\frac{r^{2}}{1+\alpha r}\right)=0.$ (26)

Figure 5:  The graph of $f_{\alpha}(r)$ for $\alpha=\frac{1}{2}$, the shape of
which is typical for $0\leq\alpha<1$. For $-1<\alpha<0$ the graph has a
vertical asymptote at $r=-\frac{1}{\alpha}$, and for $\alpha<-1$ its shape is
inverted about the $r$-axis.

The rest of the proof is similar to the proof of Theorem 5, and we only sketch
it. Note that $f_{\alpha}(0)=f^{\prime}_{\alpha}(0)=0$, and
$f^{\prime\prime}_{\alpha}(0)<0$, for all $\alpha$. Moreover (we leave the
case $\alpha=0$ to the reader),

$f_{\alpha}^{{}^{\prime\prime}}(r)=(1-\alpha)(1+\alpha
r)^{\frac{1-2\alpha}{\alpha}}-\frac{2}{(1+\alpha r)^{3}}$

stays negative for all $r>0$ when $\alpha\geq 1$. Therefore, there are no
positive solutions in this case. For $-1<\alpha<1$ the equation
$f_{\alpha}^{{}^{\prime\prime}}(r)=0$ has a single positive zero $r_{**}$ that
can be found explicitly. For $0<\alpha<1$ this means that
$f_{\alpha}^{{}^{\prime}}(r)$ decreases on $(0,r_{**})$ and strictly increases
on $(r_{**},\infty)$. Since it is negative on $(0,r_{**})$ there is a unique
$r_{*}$ on $(r_{**},\infty)$ such that $f_{\alpha}^{{}^{\prime}}(r_{*})=0$.

By a similar argument, $f_{\alpha}(r)=0$ has a unique positive root $r_{c}\geq
r_{*}\geq r_{**}>0$. The estimate $\tau_{c}\geq 2$ follows from Theorem 4, and
implies $r_{c}(1-\alpha)\geq 1$ by the first equation in (4.2). For
$-1<\alpha<0$ one should replace $\infty$ with $-\frac{1}{\alpha}$ and for
$\alpha<-1$ reverse the positive/negative and the increase/decrease in the
argument above. ∎

The case $\alpha=0$, i.e. $L(y)=\ln y$, was studied by King and Roughgarden,
who estimated $r_{c}\approx 1.793$ and $\tau_{c}\approx 2.793$. When
$\alpha\to 1$ our estimate implies that $r_{c}\to\infty$, and $\tau_{c}\to 2$,
and one can show that $r_{c}\sim-\frac{1}{\alpha}$ for large negative
$\alpha$. Graphs of $r_{c}(\alpha)$ and $\tau_{c}(\alpha)$ are depicted on
Figure 6. Biologically, this means that, all else being equal, increase in the
fitness index $\alpha$ leads to shorter times-to-go at the switch to the final
reproductive arc, and larger vegetative mass at the time of the switch.

Figure 6:  Graphs of $r_{c}(\alpha)$ (above) and $\tau_{c}(\alpha)$ (below) on
$(-\infty,1)$.

Recall that points on the switching surface that are above the cutoff curve
are the points where optimal control can switch from $u=1$ to $u=0$. In
general, these are characterized by the inequality $\textrm{Cut}>0$, where Cut
is defined in (15). We can now make this condition more explicit.

######  Lemma 4.

If $(\tau,r,y)$ is a VR junction point of an optimal trajectory then $r\leq
r_{c}$ and $\tau\leq\tau_{c}$.

###### Proof.

The $\textrm{Cut}\geq 0$ condition simplifies to

$r\leq\frac{(1+r\tau)^{\alpha}-1}{\alpha}.$

Combining this inequality with (21), we derive:

$r\tau-\frac{(1+r\tau)^{\alpha}-1}{\alpha}\left(1+\frac{1-(1+r\tau)^{-\alpha}}{\alpha}\right)\leq
0.$

The same analysis as for $f_{\alpha}(r)$ in Lemma 3 can be carried out for the
left hand side of this inequality as a function of $r\tau$. It has a single
positive zero $r_{c}\tau_{c}$, and the inequality implies that $r\tau\leq
r_{c}\tau_{c}$. By (4.2), we obtain:

$r+\frac{r^{2}}{1+\alpha r}\leq r_{c}+\frac{r_{c}^{2}}{1+\alpha r_{c}}.$

Since the left hand side is strictly increasing for positive $r$ (for
$r<-\frac{1}{\alpha}$ when $\alpha<0$) this implies $r\leq r_{c}$. The second
estimate follows from monotonicity of $\tau$ as a function of $r$ given
implicitly by (21). ∎

###  4.3 Existence of singular arcs

Note that $r_{c},\tau_{c}$ are universal functions of $\alpha$ that do not
depend on the initial values, the season length $T$, or the production
function $P(x)$. This means that the switching curve and the cutoff point do
not depend on $P(x)$ either. However, the two-dimensionalization of the
geometry described above is incomplete for general production functions. The
singular control, and therefore the singular surface, does depend on $P(x)$,
and is not a cylinder over a curve in the $r$-$\tau$ plane in general.
Nonetheless, the information on the active part of the switching surface,
provided by the last lemma, along with monotone behavior of the mass ratio $q$
along V and R arcs, is sufficient to infer the existence of singular arcs on
optimal trajectories.

The inconvenience of the smaller $r$ values being “above” the cutoff curve
when is removed by using $q=\frac{1}{r}$, which makes “above” indeed above,
see Figure 7. It follows from Lemma 4 that the region $q>q_{c}$,
$\tau>\tau_{c}$, is above the switching surface, which leads to the main
result of this section.

Figure 7:  The switching surface (blue), and the singular surface (red) for
$P(x)=x^{0.7}$ and $L(y)=2y^{0.5}$.

######  Theorem 6.

Let $L(y)=\frac{y^{\alpha}}{\alpha}$, $\alpha<1$, and $P(x)$ be concave down.
Set $t_{*}:=T-\frac{\tau_{c}}{P^{\prime}(x_{0})}$, and suppose that

$q_{c}-P^{\prime}(x_{0})t_{*}<q_{0}<q_{c}\,e^{P^{\prime}(x_{0})t_{*}}\,.$ (27)

Then, any optimal trajectory with $q_{0}$ as the initial value has a singular
arc.

###### Proof.

If the initial arc is M, we are done. If it is R then $u=0$, and
$\dot{q}=\frac{P(x_{0})}{\xi}=P^{\prime}(x_{0})$. By the lower bound in (27),
we have that $q(t_{*})>q_{c}$ and
$\tau(t_{*})=P^{\prime}(x_{0})(T-t_{*})=\tau_{c}$. This means that just before
$t_{*}$ the arc is above the switching curve, and hence is suboptimal past
$t_{*}$ by Lemma 1. Since RV switches are impossible by Theorem 3, it must
switch to an M arc before $t_{*}$.

If the initial arc is V, we can estimate:

$\dot{\xi}=\left(\frac{P(x)}{P^{\prime}(x)}\right)^{\prime}\dot{x}=\left(1-\frac{P^{\prime\prime}(x)}{P^{\prime}(x)}\frac{P(x)}{P^{\prime}(x)}\right)P(x)=\left(P^{\prime}(x)-P^{\prime\prime}(x)\xi\right)\xi\geq
P^{\prime}(x)\xi.$

Therefore, $\xi(t)\geq\xi_{0}e^{P^{\prime}(x)t}$ and
$q(t)=\frac{y_{0}}{\xi(t)}\leq q_{0}e^{-P^{\prime}(x_{0})t}$. This means that
$\tau(t)<\tau_{c}$ for $t<t_{*}$, and, by the upper bound in (27),
$q(t)<q_{c}$ for $t\geq t_{*}$. Therefore, by Lemma 4, this V arc can not
intersect the switching surface at a point where the VR switch is optimal.
Since the final arc must be R, it must switch to an M arc before $t_{*}$. ∎

######  Corollary 2.

Under the conditions of Theorem 6, if $T>\frac{\tau_{c}}{P^{\prime}(x_{0})}$
there exist optimal trajectories that include singular arcs, and if
$T>\frac{\tau_{c}+q_{c}}{P^{\prime}(x_{0})}$ then pure R trajectories are
suboptimal, and all optimal trajectories with small enough $q_{0}$ are RMR.

In other words, singular arcs (mixed growth) are a general phenomenon for
conservative ($\alpha<1$) fitness indices. On the other hand, for $\alpha\geq
1$, that corresponds to “gambler’s” strategies, we always get bang-bang
behavior.

##  5 Linear production and the mass ratio diagram

In this section we give a complete solution of the optimal allocation problem
for linear production functions and power utilities. Note that for linear
$P(x)$ we have $\xi=x$ and $q=\frac{y}{x}$, which is merely the ratio of the
reproductive to the vegetative mass of the plant. The switching and the
singular surfaces are both cylindrical in $q$-$t$ coordinates, so the
solutions can be depicted on a a two-dimensional diagram, which we call the
mass ratio diagram, see Figure 9. For $L(y)=\ln y$, this diagram was
originally described (with only partial proof) in [28]. It is somewhat
remarkable that essentially the same diagram is valid for all
$L(y)=\frac{y^{\alpha}}{\alpha}$ with $\alpha<1$, and even that the picture is
two-dimensional for $\alpha\neq 0$, as there is no obvious way to produce a
payoff functional equivalent to (3) that depends only on the mass ratio. For
$\alpha=0$ it can be done by using the additive property of logarithms.

For definiteness, we take $P(x)=x$, but one can easily reformulate the results
for $P(x)=kx$ with $k>0$ by rescaling time. The equation of the switching
curve in $q$-$t$ coordinates follows directly from (21):

$T-t=(1+\alpha+q)\frac{1-(1+\frac{T-t}{q})^{\alpha}}{\alpha}\,.$ (28)

The state equations (1) can be reduced to a single equation for the mass ratio
$q$:

$\dot{q}=1-(1+q)u.$ (29)

The singular control (7) also depends only on the mass ratio:

$u=1-\frac{1}{1-\alpha}\frac{y}{x}=1-\frac{q}{1-\alpha}.$ (30)

Figure 8:  Switching (blue) and singular (red) curves for $\alpha=-2$ in
$q$-$\tau$ coordinates, the shapes of which are typical for $\alpha<-1$. For
$-1\leq\alpha<1$ the horizontal switching asymptote is $q=0$, and for
$0\leq\alpha<1$ the horizontal singular asymptote is also $q=0$. The switching
curve is universal for all production functions, the singular curve is for
$P(x)=x$.

Integrating (29) with the singular control, we find for $\alpha\neq 0$
that333For $\alpha=0$ this simplifies to $q(t)=\frac{1}{C-t}$.

$q(t)=\frac{\alpha}{e^{\frac{\alpha}{1-\alpha}(C-t)}-1}\,.$ (31)

The constant of integration $C$ can be found from the condition that the
singular curve passes through the cutoff point, which gives
$C=T-\tau_{si}(\alpha)$, where

$\displaystyle\tau_{si}(\alpha)$
$\displaystyle=\tau_{c}(\alpha)-(1-\alpha)\ln(1+r_{c}(\alpha)\tau_{c}(\alpha))$
(32)
$\displaystyle=1+\frac{r_{c}(\alpha)}{1+\alpha\,r_{c}(\alpha)}-\frac{1-\alpha}{\alpha}\ln\Big{(}1+\alpha\,r_{c}(\alpha)\Big{)},$
(33)

and we explicitly indicated the dependence of $r_{c},\tau_{c}$ from Lemma 3 on
$\alpha$. Clearly, $\tau_{si}(\alpha)\leq\tau_{c}(\alpha)$ for $\alpha\leq 1$,
and $\tau_{si}(0)=1$. One can show that $\tau_{si}(\alpha)\geq 1$ for
$0\leq\alpha<1$. For $\alpha=-2$ we find that
$\tau_{si}=2+\sqrt{2}-3\ln(\sqrt{2}+1)\approx 0.77$, see Figure 8.

Figure 9:  The mass ratio diagram for $L(y)=\frac{y^{\alpha}}{\alpha}$ with
$\alpha<1$, and $P(x)=x$. The inactive parts of the switching and singular
curves are dashed, and the vertical lines are their asymptotes. V arcs fall
exponentially, R arcs rise with slope $1$, and M arcs follow the red curve.
Optimal trajectories starting in the labeled regions have the indicated
control structure, and MR trajectories must start on the solid part of the red
curve.

The resulting configuration of the switching and singular curves, along with
the arrangement of optimal trajectories it implies, is depicted on Figure 9.
We will now show how to determine the optimal schedules explicitly from the
initial values based on it. The proof is a more precise version of the sketch
given for the general case at the end of Section 3.3. Optimal solutions are
found analytically by integrating the state equations. For simplicity, we
state the theorem for the case of $T>q_{c}+\tau_{c}$, but the modifications
needed for simpler diagrams, where pure R trajectories occur, are
straightforward. We also leave out the modifications necessary for
$\alpha=0,-1$.

######  Theorem 7.

Let $P(x)=x,L(y)=\frac{y^{\alpha}}{\alpha}$, $\alpha<1$,
$q_{0}=\frac{y_{0}}{x_{0}}$ be the mass ratio of the initial values, and
suppose that $T>q_{c}+\tau_{c}$. Then:

(i) If
$q_{0}<\frac{\alpha}{\exp(\frac{\alpha}{1-\alpha}\left(T-\tau_{si})\right)-1}$,
then the optimal trajectory is RMR. The second switching time is
$t_{2}=T-\tau_{c}$, and the first is $t_{1}=\frac{y_{1}-y_{0}}{x_{0}}$, where
$y_{1}$ can be found by solving for $r_{1}:=\frac{x_{0}}{y_{1}}$ in

$\ln\frac{1+\alpha r_{1}}{1+\alpha
r_{0}}=\frac{\alpha}{1-\alpha}\left(T-\tau_{c}-\ln\frac{r_{1}}{r_{0}}\right);$
(34)

(ii) If
$q_{0}=\frac{\alpha}{\exp(\frac{\alpha}{1-\alpha}\left(T-\tau_{si})\right)-1}$,
then the optimal trajectory is MR. The switching time is $T-\tau_{c}$;

(iii) If
$\frac{\alpha}{\exp(\frac{\alpha}{1-\alpha}\left(T-\tau_{si})\right)-1}<q_{0}<\frac{1}{r_{c}}e^{T-\tau_{c}}$,
then the optimal trajectory is VMR. The second switching time is
$t_{2}=T-\tau_{c}$, and the first is
$t_{1}=\int\limits_{x_{0}}^{x_{1}}\frac{dx}{P(x)}$, where $x_{1}$ can be found
by solving for $r_{1}:=\frac{x_{1}}{y_{0}}$ in (34);

(iv) If $q_{0}\geq\frac{1}{r_{c}}e^{T-\tau_{c}}$, then the optimal trajectory
is VR. The switching time $t_{s}$ can be found as
$t_{s}=\ln\frac{x_{s}}{x_{0}}$ after solving for $x_{s}$ in:

$T-\ln\frac{x_{s}}{x_{0}}=\left(1+\alpha+\frac{x_{s}}{y_{0}}\right)\frac{1}{\alpha}\left(1-\left(1+\frac{x_{s}}{y_{0}}\ln\frac{x_{s}}{x_{0}}\right)^{-\alpha}\right).$
(35)

###### Proof.

One easily checks that the conditions of Theorem 3 and Lemma 1 are met in this
case.

(i) The inequality means that the trajectory starts under the singular curve.
The initial arc can not be V because V arcs drop exponentially, and no switch
to any other arc is possible then, while having a pure V trajectory
contradicts Theorem 3(i). Since $q_{0}$ is not on the singular curve the
initial arc must be R. It can either stay R to the end or switch to an M arc
only on the cutoff curve. The former possibility is suboptimal by Lemma 1, so
the second arc is M, which must switch to the final R arc by Theorem 3.
Equation (34) follows from (31) rewritten for $r=\frac{1}{q}$.

(ii) Initial V and R arcs are ruled out by the reasoning in (i) and (iii),
respectively. Hence, MR is the only possibility.

(iii) The trajectory starts above the singular curve, and hence above the
switching curve. If the initial arc were R it would rise, and hence it could
not switch to M. Therefore it must be R all the way, as RV switches are ruled
out by Theorem 3. But this is also ruled out by Lemma 1. Therefore, the
initial arc is V. It follows from the second inequality that it never crosses
the active part of the switching curve. Since it can not remain V all the way
it must switch to M when crossing the singular curve, and then to the final R
arc.

(iv) Initial R arc is ruled out as in (iii), and the initial arc is V. The
inequality implies that it never crosses the active part of the singular
surface (except at the cutoff point in the case of equality). Therefore, it
can only be VR. The equations specialize (13) and (21) to this case. ∎

##  6 Full lifetime

The model that we studied so far had the length of the safe period $T_{0}$ set
to $0$. In general, the Hamiltonian and the costate equations retain the same
form on $[T_{0},T_{0}+T]$, but on $[0,T_{0}]$ the Hamiltonian reduces to

$H(x,y,\lambda,\mu,u)=\big{(}\mu+(\lambda-\mu)u\big{)}P(x),$ (36)

and the costate equation for $\mu$ trivializes:

$\begin{cases}\dot{\lambda}=-\big{(}\mu+(\lambda-\mu)u\big{)}P^{\prime}(x)\\\
\dot{\mu}=0.\\\ \end{cases}$ (37)

This simplifies the structure of optimal trajectories on $[0,T_{0}]$. We
continue to use the labeling convention of Definition 1. The next theorem
complements Theorem 3 for the full lifetime model.

######  Theorem 8.

Let $P(x)$ and $L(y)$ be twice differentiable and $P(x)$ be positive for
$x,y>0$. Then, we have for the optimal trajectories of of problem (1), (2):

(i) Optimal trajectories have no M arcs on $[0,T_{0}]$.

(ii) If an optimal trajectory arrives at $T_{0}$ along a V arc then this arc
is initial.

(iii) If an optimal trajectory arrives at $T_{0}$ along an R arc then either
this arc is initial, or there is a single switch on $[0,T_{0}]$ from the
initial V arc.

(iv) If $P(x)$ and $L(y)$ meet all the conditions of Theorem 3 (iv) then any
optimal trajectory has no more than three switches, and one of the following
switching structures: R, VR, VMR, RMR, VRMR.

###### Proof.

(i) Along a singular arc $\dot{\lambda}=-\mu P^{\prime}(x)=\dot{\mu}=0$. Since
$P^{\prime}(x)>0$ this implies $\lambda=\mu=0$. But $\lambda$ and $\mu$ are
non-increasing over the entire interval $[0,T_{0}+T]$, so the singular arc
would have to extend all the way to $T_{0}+T$, contradicting Theorem 3 (i)
(the final arc is R).

(a) (b)

Figure 10:  Marginal values on $[0,T_{0}]$ when (a)
$\lambda(T_{0})>\mu(T_{0})$ (V arc); (b) $\lambda(T_{0})<\mu(T_{0})$ for two
different cases: VR (upper curve) and R (lower curve).

(ii) and (iii). By the above argument, $\mu(T_{0})>0$, so $\mu$ is a positive
constant on $[0,T_{0}]$. Since $(\lambda-\mu)u\geq 0$ for any admissible
control $u$, $\lambda$ is strictly monotone decreasing there by (37). If the
arc of arrival is V, then $\lambda(T_{0})\geq\mu(T_{0})$ and
$\lambda(t)>\mu(t)$ for all $t<T_{0}$. Hence, there is no switch, see Figure
10 (a). If the arc of arrival is R then $\lambda(t)=\mu(t)$ at at most one
point, see Figure 10 (b).

(iv) If the first arc on $[T_{0},T_{0}+T]$ is M or V then we have a V arc on
$[0,T_{0}]$ by (ii). Then, MR, VR, and VMR trajectories from Theorem 3 turn
into VMR, VR, and VMR, respectively. If the first arc is R then by (iii) they
either remain R on $[0,T_{0}]$, or are complemented by an initial V arc. This
means that R trajectories from Theorem 3 turn into R or VR, while RMR turn
into RMR or VRMR. ∎

The main effect of adding the safe period is, therefore, that initial V arcs
are added, or already occurring ones extended, and there is, potentially, an
additional switching time $t_{m}$ on $(0,T_{0}]$. The points $(t_{m},x_{m},y)$
form an additional switching surface that connects to the singular surface at
$t=T_{0}$, see Figure 11. If one can solve the problem on $[T_{0},T_{0}+T]$ it
is now easy to extend the solution to the entire interval.

(a) (b)

Figure 11:  Switching surface (blue), singular surface (red), and safe
switching surface (green) for $P(x)=x^{0.7}$ and $L(y)=2y^{0.5}$, shown from
two different angles.

To find the “safe” switching time $t_{m}$, note that $x$ is constant over an R
arc, and $u=0$, so $\dot{\lambda}=-\mu
P^{\prime}(x)=-\mu(T_{0})P^{\prime}(x_{m})$, where $x_{m}:=x(t_{m})$.
Therefore, on $[t_{m},T_{0}]$,

$\lambda(t)=\lambda(T_{0})+\mu(T_{0})P^{\prime}(x_{m})(T-t).$

Since $\lambda(t)$ equals $\mu(T_{0})$ at $t=t_{m}$ we derive:

$t_{m}=T_{0}-\frac{1}{P^{\prime}(x_{m})}\left(1-\frac{\lambda(T_{0})}{\mu(T_{0})}\right).$
(38)

Together with $t_{m}=\int\limits_{x_{0}}^{x_{m}}\frac{dx}{P(x)}$, which
integrates the state equation over the preceding V arc, this determines
$t_{m}$ and $x_{m}$. The time $t_{m}$ is the time when the plant is old enough
to start reproducing, and is called the age of maturity. We can show that if
the initial mass of the plant is high enough and the safe period is long
enough, the age of maturity is always positive, even if the plant has no
reproductive mass initially (as is usually the case). If the safe period is
sufficiently long, it is better to invest into growth first, to produce more
offspring later.

######  Corollary 3.

If $T_{0}P^{\prime}(x(T_{0}))\geq 1$, then the initial arc of an optimal
trajectory is always V. In particular, if $P(x)=x$, and $T_{0}>1$ then the
initial arc is V for any initial values.

###### Proof.

Since the marginal values are positive, $x(T_{0})=x_{m}$, and the growth is
purely reproductive after $t_{m}$, we must have
$t_{m}>T_{0}-\frac{1}{P^{\prime}(x_{m})}=T_{0}-\frac{1}{P^{\prime}(x(T_{0}))}>0$
by (38). ∎

The value of $x(T_{0})$ is typically easy to estimate from above in terms of
$x_{0}$ and $T_{0}$ by integrating $\dot{x}=P(x)$ on $[0,T_{0}]$. When
$P(x)=kx^{\beta}$ we have explicitly
$x(T_{0})=\left(x_{0}^{1-\beta}+(1-\beta)kT_{0}\right)^{\frac{1}{1-\beta}}$.

As before, when $P(x)=x$ and $L(y)=\frac{y^{\alpha}}{\alpha}$, the behavior of
optimal trajectories can be represented by a two-dimensional mass ratio
diagram for $q=\frac{y}{x}$, Figure 12. The maturity surface is represented by
a curve that starts from the $t$ axis and intersects the $q$ axis at the same
point as the singular curve
$q_{\textrm{sing}}=\frac{\alpha}{\exp(\frac{\alpha}{1-\alpha}\left(T-\tau_{si})\right)}$.
We denote its $t$-intercept $T_{0}-\tau_{R}$.

Figure 12:  The mass ratio diagram for the extended time model with $P(x)=x$
and $L(y)=\frac{y^{\alpha}}{\alpha}$ with $\alpha<1$. The inactive parts of
the switching and singular curves are dashed, the vertical lines are their
asymptotes. V arcs fall exponentially, R arcs rise with slope $1$, M arcs
follow the red curve. Optimal trajectories starting in the labeled regions
have the indicated control structure.

Finally, we state an analog of Theorem 7 for the case of $T>q_{c}+\tau_{c}$
and $T_{0}>\tau_{R}$, for which the proof is analogous. The modifications for
shorter times, when pure R and RMR trajectories also occur, are
straightforward.

######  Theorem 9.

Let $P(x)=x,L(y)=\frac{y^{\alpha}}{\alpha}$, and $q_{0}=\frac{y_{0}}{x_{0}}$
be the ratio of initial values. Then:

(i) If $q_{0}<\frac{\alpha
e^{T_{0}}}{\exp(\frac{\alpha}{1-\alpha}\left(T-\tau_{si})\right)-1}$, then the
optimal trajectory is VRMR.

(ii) If $\frac{\alpha
e^{T_{0}}}{\exp(\frac{\alpha}{1-\alpha}\left(T-\tau_{si})\right)-1}\leq
q_{0}<\frac{1}{r_{c}}e^{T_{0}+T-\tau_{c}}$, then the optimal trajectory is
VMR.

(iii) If $q_{0}\geq\frac{1}{r_{c}}e^{T_{0}+T-\tau_{c}}$, then the optimal
trajectory is VR.

Acknowledgments: This work was conceived during the summer 2019 REU program at
the University of Houston-Downtown, and is funded by the National Science
Foundation grant 1560401. The authors are grateful to the anonymous reviewers
for clarifying the relation between density dependence and fitness and many
other helpful suggestions that greatly improved the paper.

## References

  * [1] Amir, S., Cohen, D.: Optimal reproductive efforts and the timing of reproduction of annual plants in randomly varying environments, Journal of Theoretical Biology 147 (1990) 17–42.
  * [2] Argasinski, K., Broom, M., The nest site lottery: how selectively neutral density dependent growth suppression induces frequency dependent selection, Theoretical Population Biology, 90 (2013) 82–90.
  * [3] Argasinski, K., Rudnicki, R., Nest site lottery revisited: towards a mechanistic model of population growth suppressed by the availability of nest sites, Journal of Theoretical Biology, 420 (2017) 279–289.
  * [4] Argasinski, K., Rudnicki, R., From nest site lottery to host lottery: continuous model of growth suppression driven by the availability of nest sites for newborns or hosts for parasites and its impact on the selection of life history strategies, Theory in Biosciences, 139 (2020) 171–188.
  * [5] Cohen, D., Optimizing reproduction in a randomly varying environment, Journal of Theoretical Biology, 12 (1966) 19–129.
  * [6] Cohen, D., Maximizing final yield when growth is limited by time or by limiting resources, Journal of Theoretical Biology, 33 (1971) no. 2, 299–307.
  * [7] Brown, J., Venable, D., Life history evolution of seed-bank annuals in response to seed predation, Evolutionary Ecology, 5 (1991) 12–29.
  * [8] Bell, G., Lechowicz, M., Schoen, D., The Ecology and Genetics of Fitness in Forest Plants III, Journal of Ecology , 79 (1991) no. 3 pp. 697–713.
  * [9] Clark, C., Poulsen, J., Levey, D., Osenberg, C., Are plant populations seed limited? A critique and meta-analysis of seed addition experiments, The American Naturalist, 170 (2007) no. 1, 128–142.
  * [10] De Lara, M., Mum, why do you keep on growing? Impacts of environmental variability on optimal growth and reproduction allocation strategies of annual plants, Journal of Mathematical Biology, 52 (5) (2006) 633–666.
  * [11] Dempster, E., 1955. Maintenance of genetic heterogeneity, Quantitative Biology, 20 (1955) 2–32.
  * [12] Denholm, J., Necessary condition for maximum yield in a senescing two-phase plant, Journal of Theoretical Biology, 52(1) (1975) no. 2, 251–254.
  * [13] Ejsmond, A., Kozlowski, J., Ejsmond, M., Probing of mortality rate by staying alive: The growth-reproduction trade-off in a spatially heterogeneous environment, Functional Ecology, 33 (2019) 2327–2337.
  * [14] Engen, S., Saether, B. (2017) $r$-and $K$-selection in fluctuating populations is determined by the evolutionary trade-off between two fitness measures: growth rate and lifetime reproductive success, Evolution, 71(1) (2017) 167–173.
  * [15] Fox, G., Evolutionary Ecology, Annual Review of Ecology and Systematics, 24, (1992) 482–499.
  * [16] Frank, S., Natural selection I. Variable environments and uncertain returns on investment, Journal of Evolutionary Biology, 24 (2011) 2299–2309.
  * [17] Frank, S., Slatkin, M., Evolution in a variable environment, The American Naturalist, 136 (1990) 244–260.
  * [18] Gillespie, J., Natural selection for variances in offspring numbers: a new evolutionary principle, The American Naturalist, 111 (1977) 1010–1014.
  * [19] Goudriaan, J., A Family of Saturation Type Curves, Especially in Relation to Photosynthesis, Annals of Botany, 43 (1979) 783–785.
  * [20] Harper, J., Population biology of plants, Academic Press, London, 1977.
  * [21] Hui, C., Carrying capacity, population equilibrium, and environment’s maximal load, Ecological Modelling, 192 (2006) no. 1–2, 317–320.
  * [22] Hulme, P., Benkman, C., Granivory, In: Plant-animal interactions: an evolutionary approach, Herrera, C., Pellmyr, O. (eds), Blackwell Science, NY, 2002, 132–154.
  * [23] Inouye, R., Density-dependent germination response by seeds of desert annuals, Oecologia, 46 (1980) 235–238.
  * [24] Ioslovich, I., Gutman, P.-O., On the botanic model of plant growth with intermediate vegetative-reproductive stage, Theoretical Population Biology, 68 (2005) 147–156.
  * [25] Iwasa Y., Roughgarden, J., Shoot/root balance of plants: Optimal growth of a system with many vegetative organs, Theoretical Population Biology, 25 (1984) no. 1, 78–105.
  * [26] Johansson, J., Brännström, A., Metz, J., Dieckmann, U., Twelve fundamental life histories evolving through allocation-dependent fecundity and survival, Ecology and Evolution, 8(6) (2018) 3172–3186.
  * [27] Kelly, D., Sork, V., Mast seeding in perennial plants: why, how, where? Annual Review of Ecology and Systematics, 33 (2002) pp. 427–447.
  * [28] King, D., Roughgarden, J., Graded allocation between vegetative and reproductive growth for annual plants in growing seasons of random length, Theoretical Population Biology, 22 (1982) no. 1, 1–16.
  * [29] Lande, R., Engen, S., Saether, B., An evolutionary maximum principle for density-dependent population dynamics in a fluctuating environment, Philosophical Transactions of the Royal Society B, 364 (2009) 1511–1518.
  * [30] Lenhart, S., Workman, J., Optimal control applied to biological models, Chapman & Hall/CRC, Boca Raton, FL, 2007.
  * [31] Lewontin, R., Cohen, D., On population growth in a randomly varying environment, PNAS, 62(4) (1969) 1056–1060.
  * [32] Linhart, Y., Density-Dependent Seed Germination Strategies in Colonizing Versus Non-Colonizing Plant Species, Journal of Ecology, 64 (1976) no. 1, pp. 375–380.
  * [33] MacArthur, R., Wilson, E., The theory of island biogeography. Princeton University Press, Princeton, 1967.
  * [34] Mironchenko, A., Kozlowski, J., Optimal allocation patterns and optimal seed mass of a perennial plant, Journal of Theoretical Biology, 354 (2014) 12–24.
  * [35] Ooi, R., Whelan, J., Auld, T., Persistence of obligate-seeding species at the population scale: effects of fire intensity, fire patchiness and long fire-free intervals, International Journal of Wildland Fire, 15 (2006) pp. 261–269.
  * [36] Perrin, N., Sibly, R., Dynamic Models of Energy Allocation and Investment, Annual Review of Ecology and Systematics, 24, (1993) 379–410.
  * [37] Pianka, E., On $r$\- and $K$-selection, The American Naturalist, 104 (1970) no. 940, 592–597.
  * [38] Post, D., Palkovacs, E., Eco-evolutionary feedbacks in community and ecosystem ecology: interactions between the ecological theatre and the evolutionary play, Philosophical Transactions of the Royal Society B, 364, (2009) 1629–1640.
  * [39] Reproductive Allocation in Plants, Edward G. Reekie and Fakhri A. Bazzaz (eds.), Academic Press, Burlington, MA, 2005.
  * [40] Roff, D., Defining fitness in evolutionary models, Journal of Genetics, 87 (2008) 339–348.
  * [41] Seger J., Brockmann H., What is bet-hedging? Oxford surveys in evolutionary, 4 (1987) pp. 182–211.
  * [42] Seierstad, A., Sydsaeter, K., Optimal control theory with economic applications, Elsevier, New York, 1986.
  * [43] Slatkin, M., Hedging one’s evolutionary bets, Nature, 250 (1974) 704–705.
  * [44] Takahashi,D., Yamauchi, A., Optimal defense schedule of annual plants against seasonal herbivores, American Naturalist 175(5) (2010) 538–550.
  * [45] Venable, D., Bet hedging in a guild of desert annuals, Ecology, 88(5) (2007) pp. 1086–1090.
  * [46] Vincent, T., Pulliam, H., Evolution of life history strategies for an asexual annual plant model, Theoretical Population Biology, 17 (1980), no. 2, 215–231.
  * [47] Wenk, E., Falster, D., Quantifying and understanding reproductive allocation schedules in plants, Ecology and Evolution, 5 (2015) no. 23, 5521–5538.
  * [48] Zhang, H., Cheng, J., Xiao, Z. et al., Effects of seed abundance on seed scatter-hoarding of Edward’s rat (Leopoldamys edwardsi Muridae) at the individual level, Oecologia, 158 (2008) 57–63.
  * [49] Ziolko, M., Kozlowski,J., Some optimization models of growth in biology, IEEE Transactions Automatic Control, 40(10) (1995) 1779–1783.