
The fitness of an annual plant can be thought of as how much fruit is produced by the end of its growing season. Working under the assumption that annual plants grow to maximize fitness, we use optimal control theory to understand this process. We introduce a model for resource allocation in annual plants that extends classical work by Iwasa and Roughgarden to a case where both carbohydrates and mineral nutrients are allocated to shoots, roots, and fruits. We use optimal control theory to determine the optimal resource allocation strategy for the plant throughout its growing season as well as develop a numerical scheme to implement the model. We find that fitness is maximized when the plant undergoes a period of mixed vegetative and reproductive growth prior to switching to reproductive-only growth at the end of the growing season. Our results further suggest that what is optimal for an individual plant is highly dependent on initial conditions, and optimal growth has the effect of driving a wide range of initial conditions toward common configurations of biomass by the end of a growing season.
Citation: David McMorris, Glenn Ledder. Optimal allocation of two resources in annual plants[J]. Mathematical Biosciences and Engineering, 2025, 22(6): 1464-1516. doi: 10.3934/mbe.2025055
[1] | Agustín Gabriel Yabo, Jean-Baptiste Caillau, Jean-Luc Gouzé . Optimal bacterial resource allocation: metabolite production in continuous bioreactors. Mathematical Biosciences and Engineering, 2020, 17(6): 7074-7100. doi: 10.3934/mbe.2020364 |
[2] | Shahab Shamshirband, Javad Hassannataj Joloudari, Sahar Khanjani Shirkharkolaie, Sanaz Mojrian, Fatemeh Rahmani, Seyedakbar Mostafavi, Zulkefli Mansor . Game theory and evolutionary optimization approaches applied to resource allocation problems in computing environments: A survey. Mathematical Biosciences and Engineering, 2021, 18(6): 9190-9232. doi: 10.3934/mbe.2021453 |
[3] | Zehan Liu, Daoxin Qiu, Shengqiang Liu . A two-group epidemic model with heterogeneity in cognitive effects. Mathematical Biosciences and Engineering, 2025, 22(5): 1109-1139. doi: 10.3934/mbe.2025040 |
[4] | Xuyin Wang, Weiguo Liu, Lu Li, Peizhen Zhao, Ruifeng Zhang . Resource dependent scheduling with truncated learning effects. Mathematical Biosciences and Engineering, 2022, 19(6): 5957-5967. doi: 10.3934/mbe.2022278 |
[5] | Yanpei Liu, Yunjing Zhu, Yanru Bin, Ningning Chen . Resources allocation optimization algorithm based on the comprehensive utility in edge computing applications. Mathematical Biosciences and Engineering, 2022, 19(9): 9147-9167. doi: 10.3934/mbe.2022425 |
[6] | Xinru Zhou, Xinmiao Rong, Meng Fan, Josué-Antonio Nescolarde-Selvaa . Stoichiometric modeling of aboveground-belowground interaction of herbaceous plant. Mathematical Biosciences and Engineering, 2019, 16(1): 25-55. doi: 10.3934/mbe.2019002 |
[7] | Semu Mitiku Kassa . Three-level global resource allocation model for HIV control: A hierarchical decision system approach. Mathematical Biosciences and Engineering, 2018, 15(1): 255-273. doi: 10.3934/mbe.2018011 |
[8] | Yafei Li, Yuxi Liu . Multi-airport system flight slot optimization method based on absolute fairness. Mathematical Biosciences and Engineering, 2023, 20(10): 17919-17948. doi: 10.3934/mbe.2023797 |
[9] | Ali Gharouni, Lin Wang . Modeling the spread of bed bug infestation and optimal resource allocation for disinfestation. Mathematical Biosciences and Engineering, 2016, 13(5): 969-980. doi: 10.3934/mbe.2016025 |
[10] | Yu Shen, Hecheng Li . A new differential evolution using a bilevel optimization model for solving generalized multi-point dynamic aggregation problems. Mathematical Biosciences and Engineering, 2023, 20(8): 13754-13776. doi: 10.3934/mbe.2023612 |
The fitness of an annual plant can be thought of as how much fruit is produced by the end of its growing season. Working under the assumption that annual plants grow to maximize fitness, we use optimal control theory to understand this process. We introduce a model for resource allocation in annual plants that extends classical work by Iwasa and Roughgarden to a case where both carbohydrates and mineral nutrients are allocated to shoots, roots, and fruits. We use optimal control theory to determine the optimal resource allocation strategy for the plant throughout its growing season as well as develop a numerical scheme to implement the model. We find that fitness is maximized when the plant undergoes a period of mixed vegetative and reproductive growth prior to switching to reproductive-only growth at the end of the growing season. Our results further suggest that what is optimal for an individual plant is highly dependent on initial conditions, and optimal growth has the effect of driving a wide range of initial conditions toward common configurations of biomass by the end of a growing season.
Plant life history theory is generally concerned with the strategies that plants employ for survival and reproduction, as well as how these processes influence population dynamics. The question of how plants allocate resources, e.g., C, N, P, and in some models, biomass, and what drives these allocation rules, is central to this pursuit. There are several schools of thought that seek to provide a framework for understanding these allocation patterns. One school of thought suggests that allocation rules should ultimately be the result of natural selection, and so resource allocation should optimize fitness in some sense, for example by maximizing fruits [1] or overall biomass [2]. Another framework views biomass allocation as following certain allometric scaling relationships [3], which has been shown to be compatible with the optimization framework [4]. Yet another views allocation not through the lens of an individual organism or specific genome but rather from a game-theoretical perspective in which allocation rules are driven by competition and follow an evolutionary stable strategy [5]. It has also been suggested that allocation may be an emergent property of the many processes involved in plant growth rather than the result of a central allocation process [6]. For a more complete review of allocation theory, we refer the reader to [7] or [8].
In this paper, we take the viewpoint that resource allocation, in annual plants specifically, should serve to optimize overall fitness. It is important to note that, while patterns of growth consistent with optimal allocation have been observed to some extent [4], even if this is not universally true, it is still important to have a theory of optimal growth for comparison with the observed behavior. Whereas previous work has focused primarily on optimal allocation of a single resource, be it carbon or biomass, and ignored the role of mineral nutrients, our work seeks to develop a theory that acknowledges the importance of both carbon and mineral nutrients. This is in line with the functional equilibrium hypothesis, which states that optimal growth occurs when resources are allocated in such a way that no single resource is any more limiting than any other. This idea has its origin in economics [9] and has been studied across a wide range of experiments [10], although allometry may also play a role [11].
This work is primarily an extension of the classical work by Iwasa and Roughgarden [1], who considered a model for photosynthate (C) allocation in annual plants. They assumed that photosynthate production depends on both shoots and roots, and that allocation is split between shoots, roots, and fruits so as to maximize fruit yield. They used optimal control theory to determine that fruit yield is maximized by a three-phase growth pattern, characterized by an initial phase of shoot-only or root-only growth, a period of balanced growth during which shoots and roots grow simultaneously, and ultimately a period of fruit-only growth at the end of the growing season.
We present a model that extends the work of Iwasa and Roughgarden to a case that directly incorporates the role of mineral nutrients. In particular, we model the allocation of carbohydrates and mineral nutrients in an annual plant with the objective of optimizing fruit yield and use optimal control theory to determine the allocation strategies. In this paper, we will consider a case where shoots and roots require both resources, and, for simplicity, fruits only require carbohydrates, as in [1]. This simplification allows us to obtain mathematical results that provide an important proof of concept for the viability of this type of model, which we intend to use in the future to guide more general two-resource models.
We find, among other things, that this addition results in four phases of growth, rather than three as seen in [1]. This additional phase consists of a period of mixed vegetative/reproductive growth, during which the fruits and shoots grow simultaneously. There has been some criticism of the use of optimal control models for carbon allocation because they typically predict a bang-bang control with a single switch between vegetative and reproductive growth, whereas this separation of function is not biologically realistic [12]. We find the emergence of this additional phase in our model to be evidence for its validity, even with the simplifying assumptions we are making. Furthermore, our results indicate that what is optimal for one plant may not be optimal for another, and optimal growth is largely dependent on initial conditions. We also conclude that this range of optimal strategies may have the effect of reducing population-level variation, thus driving a population toward common sizes and optimal yields. The question remains, however, as to whether plants actually have this degree of strategic plasticity.
We will use the notation x(t)=(x1(t),x2(t),…,xn(t)) for the function x:R→Rn and g(t,x,u)=(g1(t,x,u),g2(t,x,u),…,gn(t,x,u)) for the function g:R×Rn×Rm→Rn. Of particular importance to us are a well-known class of optimal control problems of the form
maxu∫t1t0f(t,x,u)dtsubject to:x′(t)=g(t,x,u),x(t0)=x0subject to:ak≤uk≤bkfor k=1,…,m. | (2.1) |
The solution process is formalized in Pontryagin's maximum principle [13,14,]. Here, we form a Hamiltonian with n adjoints λi, one for each of the states xi:
H(t,x,u,λ)=f(t,x,u)+λ(t)⋅g(t,x,u), | (2.2) |
where each adjoint satisfies the following:
λ′i(t)=−∂H∂xi(t,x∗,u∗,λ),λi(t1)=0for i=1,…,n. | (2.3) |
An optimal pair (x∗,u∗) is characterized by a set of necessary conditions for optimality (2.4). We refer the reader to [13] for a thorough treatment of this class of optimal control problems.
{uk=ak,if ∂H∂uk<0ak≤uk≤bk,if ∂H∂uk=0uk=bk,if ∂H∂uk>0 | (2.4) |
The control problem we will study is an extension of this class of problems, with additional constraints on the controls. We will discuss how the necessary conditions change when they are introduced, as well as include the derivations in Appendix A.
We now turn our attention to the model for resource allocation in annual plants we will discuss as an extension of the work of Iwasa and Roughgarden [1] to a scenario where the growth of the roots and shoots relies on two resources, rather than one. For the sake of convenience, we will use the terms carbon and nitrogen to refer to more complicated classes of carbohydrates and mineral nutrients. As a simplifying assumption, we assume that fruits only require carbon for growth.
This assumption makes the mathematical analysis more feasible, and aligns biologically with plants that encase small seeds in large carbon-rich fruits. While this may match our colloquial definition of fruits, e.g., watermelons, it is not an accurate representation of many types of annual plants for which the term fruits refers to seeds, nuts, etc., and for which nitrogen content is not negligible.
We consider an annual plant with three organs: shoots, roots, and fruits. Shoots consist of all above-ground vegetative biomass, and roots consist of all below-ground vegetative biomass. Fruits refer to any reproductive biomass, be it in the form of colloquial fruits, nuts, seeds, etc. The biomass of each organ is measured in units of carbon, which is used as a catch-all term for carbohydrates produced by the shoots. The functions S(t),R(t), and F(t) give the biomass of shoots, roots, and fruits, respectively, at time t throughout a fixed growing season [0,T]. We assume that the plant relies on two resources, carbon and nitrogen, though we use these terms loosely to refer to more complicated classes of carbohydrates and soil nutrients. Carbon is fixed by the shoots at a rate of C(S), and nitrogen is absorbed by the roots at a rate of N(R). Note that here, we made a simplifying assumption, as in reality the rate of carbon fixation depends on both shoots and roots via transpiration. Throughout the growing season, fractions of carbon uSC(t),uRC(t), and uFC(t) and fractions of nitrogen uSN(t),uRN(t), and uFN(t) are allocated to the shoots, roots, and fruits, respectively, at time t, and we assume that the plant can freely allocate resources. The resources pass through a synthesizing unit (SU) in each organ, where they are converted into biomass. We use Kooijman's parallel complementary synthesizing unit (PCSU) function from [15], employing the same simplification seen in [7]HY__HY, Appendix A], given by
g(A,B)=A2B+AB2A2+AB+B2, | (3.1) |
to provide the rate of tissue production when resources are provided at rates A and B. We can view these rates as representing the maximum rate of tissue production given the full utilization of that particular resource. Therefore, when resources are available in the stoichiometric ratio, we have that A=B. As the tissue is measured in units of carbon, we need conversion factors νS,νR, and νF, which give the fixed stoichiometric C:N ratios in the shoots, roots, and fruits, respectively. We specify initial conditions S(0)=S0,R(0)=R0, and F(0)=0, which leads to the following model:
dSdt=g(uSCC,νSuSNN),S(0)=S0 | (3.2) |
dRdt=g(uRCC,νRuRNN),R(0)=R0 | (3.3) |
dFdt=g(uFCC,νFuFNN),F(0)=0 | (3.4) |
As we are considering the case where fruits only require carbon, we implement this assumption by taking νF→∞, and note that for any nonzero uFNN, we have the limit
limνF→∞g(uFCC,νFuFNN)=uFCC. | (3.5) |
This leads to the following simplified model, which will form the basis for our analysis.
dSdt=g(uSCC,νSuSNN),S(0)=S0 | (3.6) |
dRdt=g(uRCC,νRuRNN),R(0)=R0 | (3.7) |
dFdt=uFCC,F(0)=0 | (3.8) |
Note that we suppress the arguments in most functions for convenience. This model is shown schematically in Figure 1.
Because the functions uFC,uSC,uRC,uSN, and uRN represent fractions of carbon and nitrogen we impose several restrictions. First, we require that each function be piecewise continuous and bounded between 0 and 1. Furthermore, as we assume full utilization of each resource, we assume that for all times t∈[0,T]
uFC+uSC+uRC=1=uSN+uRN. | (3.9) |
Additionally, to be biologically realistic, we assume that the plant's capacity to collect resources increases continuously with biomass, meaning that we require both C and N to be continuously differentiable, and
dCdS>0anddNdR>0. | (3.10) |
We also assume that the plant experiences possibly diminishing returns with increased biomass, meaning that
d2CdS2≤0andd2NdR2≤0. | (3.11) |
Before we introduce the optimal control framework we will use for determining the optimal resource allocation strategy, there are several important features of the PCSU that merit discussion. As introduced in (3.1), for resource fluxes A and B, the PCSU function
g(A,B)=A2B+AB2A2+AB+B2, | (3.1 revisited) |
gives the rate of tissue production. Note that when only one resource is present, the tissue production rate is zero as g(A,0)=0=g(0,B) for nonzero A and B. Additionally, this is still true when neither resource is present. In particular, writing (3.1) in polar coordinates allows us to show that g(A,B)→0 as (A,B)→(0,0). Because of this, we will take g to be the continuous extension of this function to the origin such that g(0,0)=0.
Now, by changing variables to either zB=BA or zA=AB, we can rewrite g as
g(A,B)=AG(zB)=BG(zA), | (3.12) |
where
G(z)=z(1+z)1+z+z2. | (3.13) |
The fact that we can represent the PCSU function in this manner will aid in our analysis by restricting the nonlinearities present to a function of a single variable. By the above discussion, we observe that although zB and zA may be undefined when either A or B or both are zero, both AG(zB) and BG(zA) are zero when at least one of A or B is zero. It is important to note, however, that while g is continuous at the origin, G(zB) and G(zA) are not. Additionally, it will facilitate later analysis of the model to note that
∂g∂A(A,B)=∂∂ABG(zA)=G′(zA), | (3.14) |
so long as (A,B)≠(0,0), and, recalling that g(0,0)=0,
∂g∂A(0,0)=limh→0g(0+h,0)−g(0,0)h=limh→00h=0, | (3.15) |
although neither G′(zB) nor G′(zA) is continuous at the origin. We can obtain ∂g∂B(A,B) in the same manner, and so we have
∂g∂X(A,B)={G′(zX),(A,B)≠(0,0)0,(A,B)=(0,0)for X∈{A,B} | (3.16) |
where
G′(z)=1+2z(1+z+z2)2. | (3.17) |
Finally, it is worth noting that both G and G′ are bounded between 0 and 1 and that
G(0)=0,G′(0)=1,limz→∞G(z)=1,limz→∞G′(z)=0. | (3.18) |
There are two identities related to the PCSU function and its derivatives that will be repeatedly cited in later analysis. We will start with the following notation:
G2(z)=G′(1/z)=z3(2+z)(1+z+z2)2. | (3.19) |
Employing (3.19), we have the following. We omit the proofs as each can be verified directly.
PCSU Identity 1.
G(z)−zG′(z)=G2(z). | (3.20) |
PCSU Identity 2.
G′2(z)=−zG″(z). | (3.21) |
In this section, we will translate the original problem of finding the growth trajectory to maximize fruit yield to an optimal control problem. Noting that
F(T)=∫T0dFdtdt+F(0)=∫T0uFC(t)C(S)dt, | (3.22) |
we can reframe our goal as maximizing the right-hand side of (3.22). By imposing the previously discussed constraints that the fractions of each resource, henceforth referred to as the controls, be bounded in [0,1], and fractions of the same resource sum to unity, along with the differential equations for the biomass of each vegetative organ, we can associate the optimal growth trajectory with the solution to the following optimal control problem.
maxu∫T0uFCCdtsubject to:uiC≥0, ujN≥0for i=0,1,2, j=1,2subject to:uFC+uSC+uRC=1=uSN+uRNsubject to:dSdt=g(uSCC,νSuSNN),S(0)=S0subject to:dRdt=g(uRCC,νRuRNN),R(0)=R0 | (3.23) |
Note that here, we only require that the controls be non-negative, because implicit in the combination of those constraints and the equality constraints is the requirement that each control must also be less than one.
We will solve the optimal control problem (3.23) using a set of necessary conditions that must be satisfied by the solution. The presence of the two equality constraints makes this problem non-standard, so we take the time to derive the necessary conditions for this type of problem. This derivation is found in Appendix A. Since we can essentially think of the carbon controls and the nitrogen controls separately, the necessary conditions for (3.23) are a combination of the conditions for the two types of problems discussed in Appendix A. We begin by forming a Hamiltonian with two piecewise differentiable adjoints, λS(t) and λR(t):
H=uFCC+λSg(uSCC,νSuSNN)+λRg(uRCC,νRuRNN). | (3.24) |
The necessary conditions for optimality are as follows:
{uiC=0if∂H∂uiC<∂H∂ujC, for all j≠iuiC=1if∂H∂uiC>∂H∂ujC, for all j≠i0≤uiC≤1if∂H∂uiC=∂H∂ujC, for any j≠iuiN=0,ujN=1if∂H∂uiN<∂H∂ujN0≤uSN,uRN≤1if∂H∂uSN=∂H∂uRN | (3.25) |
What this essentially says is that all the carbon is being allocated to the organ such that the partial derivative of the Hamiltonian with respect to that organ's carbon control is larger than the partial derivatives of the Hamiltonian with respect to the other two carbon controls. Likewise for the nitrogen controls. Note also that the Hamiltonian must be constant along the optimal trajectory because the optimal control problem (3.23) is autonomous [13].
As mentioned in Section 3.2, we can rewrite the PCSU function g by a change of variables. In particular, we can rewrite (3.6), (3.7), and (3.8) via (3.12) and the substitutions
zSC=uSCCνSuSNN | (3.26) |
zRC=uRCCνRuRNN. | (3.27) |
This results in the following system
dSdt=νSuSNNG(zSC),S(0)=S0 | (3.28) |
dRdt=νRuRNNG(zRC),R(0)=R0 | (3.29) |
dFdt=uFCC,F(0)=0 | (3.30) |
where G is given by (3.13). Therefore, the Hamiltonian can be rewritten as
H=uFCC+λSνSuSNNG(zSC)+λRνRuRNNG(zRC). | (3.31) |
Now, as we previously mentioned, we can essentially think of the carbon controls and nitrogen controls separately. The above formulation provides an avenue for restricting the appearance of the carbon controls to the argument of G by letting us write
g(Carbon Flux,Nitrogen Flux)=(Nitrogen Flux)⋅G(Carbon FluxNitrogen Flux). |
When working directly with the nitrogen controls, it will be useful to write
g(Carbon Flux,Nitrogen Flux)=(Carbon Flux)⋅G(Nitrogen FluxCarbon Flux). |
With this in mind, we can make the substitutions
zSN=νSuSNNuSCC | (3.32) |
zRN=νRuRNNuRCC, | (3.33) |
and rewrite (3.6), (3.7), and (3.8) as
dSdt=uSCCG(zSN),S(0)=S0 | (3.34) |
dRdt=uRCCG(zRN),R(0)=R0 | (3.35) |
dFdt=uFCC,F(0)=0. | (3.36) |
As before, the Hamiltonian can be rewritten as
H=uFCC+λSuSCCG(zSN)+λRuRCCG(zRN). | (3.37) |
Using (3.31) and (3.37), we can compute the following partial derivatives:
∂H∂uFC=C | (3.38) |
∂H∂uSC=λSνSuSNNG′(zSC)∂∂uSCzSC=λSCG′(zSC) | (3.39) |
∂H∂uRC=λRνRuRNNG′(zRC)∂∂uRCzRC=λRCG′(zRC) | (3.40) |
∂H∂uSN=λSuSCCG′(zSN)∂∂uSNzSN=λSνSNG′(zSN) | (3.41) |
∂H∂uRN=λRuRCCG′(zRN)∂∂uRNzRN=λRνRNG′(zRN). | (3.42) |
It is important to recall that by (3.16), these formulations of ∂g∂A(A,B) and ∂g∂B(A,B) are only valid so long as (A,B)≠(0,0), in which case these partial derivatives vanish.
The last piece of the optimal control framework concerns the adjoints λS and λR. By the derivation of the necessary conditions in Appendix A, we have that
λ′S=−∂H∂S,λS(T)=0 | (3.43) |
λ′R=−∂H∂R,λR(T)=0. | (3.44) |
Making use of (3.31) and (3.37), we can express (3.43) and (3.44) as follows:
λ′S=−CS[uFC+λSuSCG′(zSC)+λRuRCG′(zRC)] | (3.45) |
λ′R=−NR[λSνSuSNG′(zSN)+λRνRuRNG′(zRN)] | (3.46) |
Note that here, we use the notation CS=C′(S) and NR=N′(R) for the sake of brevity and to reserve the prime notation for cases where either the derivative is taken with respect to time or when the argument of the function is made explicit.
In this section, we will discuss the structure of the solution to (3.23). Generally speaking, the solution exhibits a four-phase structure. First, there is an initial phase of vegetative growth in which the plant addresses deficiencies in either shoots or roots. Second, the plant undergoes a phase of balanced growth when both the shoots and roots are growing. Third, there is a phase of mixed vegetative/reproductive growth during which the shoots continue to grow and the fruits begin growing simultaneously. The plant completes the growing season with a period of reproductive growth, during which only the fruits are growing.
We will begin with the analysis of the four-phase structure of the optimal solution. Because we know that the adjoints, λS and λR, vanish at the end of the growing season, we use that as the starting point for our analysis and proceed in reverse.
We begin by showing that there exists a switching time after which all carbon is allocated to fruit production. We will refer to the interval between the switching time and T the final interval. Note that by (3.43) and (3.44), we have that λS(T)=0=λR(T). Furthermore, recall that G′ is bounded. Therefore, writing C(S(T))=C∗, we have by (3.38), (3.39), and (3.40) that
∂H∂uFC(T)=C∗>0 | (4.1) |
∂H∂uSC(T)=0 | (4.2) |
∂H∂uRC(T)=0. | (4.3) |
This implies that at t=T
∂H∂uFC(t)>∂H∂uSC(T)and∂H∂uFC(t)>∂H∂uRC(T). | (4.4) |
Therefore, by (3.25), we have that uFC(T)=1 and uSC(T)=0=uRC(T). That is, at the end of the growing season, the plant is allocating all of the available carbon to the fruits, which in turn means that only the fruits are growing at this time.
Now, because both λS and λR are continuous, λS(T)=0=λR(T), and G′ is bounded, there must be some ε>0 such that for all t in [T−ε,T], we have that (4.4) still holds. So, we have the existence of a (potentially small) interval of fruit-only growth. Now, writing CS(S(T))=C∗S, during this interval we have by (3.45) that
λ′S=−C∗S | (4.5) |
because here, we have uFC=1,uSC=0=uRC and G′ is bounded. Because λS(T)=0, this implies that during this interval
λS(t)=C∗S⋅(T−t). | (4.6) |
Additionally, because during this interval we have a scenario where for i=1,2 either only uiC is zero or both uiC and uiN are zero, we have by (3.46) that
λ′R=0. | (4.7) |
This is because either uiC alone is zero, in which case ziN→∞ and by (3.18) we have that G′(ziN)→0, or in the case that both uiC and uiN are zero, then uiNG′(ziN)=0 because G′ is bounded. In either case, we get (4.7). So, because λR(T)=0, we have that the following holds throughout the interval:
λR(t)=0. | (4.8) |
Therefore, because G′ is bounded, we have by (3.42) that
∂H∂uRN=0. |
Since (3.41) is non-negative in this interval, it must be the case that
∂H∂uRN≤∂H∂uSN. | (4.9) |
If this inequality were strict, then by (3.25), we would have uSN=1. In this case, because uSC=0, we have that zSN→∞, and so by (3.18), we have that G′(zSN)→0. This, however, implies that
∂H∂uSN=0=∂H∂uRN, |
and so the inequality (4.9) cannot be strict, and both partial derivatives are zero in this interval. In particular, G′(zSN)=0, which means that zSN→∞ and so zSC=0, which by (3.18) implies that G′(zSC)=1. In summary, during this interval, we have
∂H∂uFC=C∗,∂H∂uSC=λSC∗,∂H∂uRC=∂H∂uSN=∂H∂uRN=0. | (4.10) |
Now, since λS is decreasing by (4.5), we see from (4.10) that the partial derivatives will maintain the same ordering as long as λS<1. We define the switching point t∗ to be the time when the plant switches to fruit-only growth. We can identify the switching point as the time when λS(t)=1, and we can use (4.6) to define t∗ implicitly by the equation
C∗S⋅(T−t∗)=1. | (4.11) |
By solving for t∗ in (4.11), we obtain
t∗=T−1C∗S. | (4.12) |
In doing so, note that we have extended this period of fruit-only growth from the interval [T−ε,T] to [t∗,T]. Note also that during the final interval, we have by (3.24) that
H=C, |
and because the Hamiltonian must be constant along the optimal trajectory, we have that
H=C∗ | (4.13) |
throughout the growing season.
We will now continue backwards to show that, prior to the final interval of fruit-only growth, we have a period of mixed vegetative/reproductive growth, during which both the shoots and fruits grow simultaneously. This interval will be referred to as the penultimate interval. Recall that at the beginning of the final interval, we have that λS=1, and so (4.10) becomes
∂H∂uFC=∂H∂uSC=C,∂H∂uRC=∂H∂uSN=∂H∂uRN=0. | (4.14) |
In order to apply (3.25) we need to determine the ordering of the partial derivatives of H in an interval immediately prior to t∗. We will accomplish this by a sequence of three lemmas, the proofs of which are included in Appendix B.1. We begin with the following inequality.
Lemma 1. ∂H∂uRC<∂H∂uFC in an open interval immediately prior to t∗.
Note that, regardless of where ∂H∂uSC falls in the order of partial derivatives, the strict inequality in Lemma 1 means that there is no root growth during this period of time immediately prior to t∗. Next, we have the following, which indicates that allocation to shoots and to fruits at this point in the growing season are equally important to overall fitness.
Lemma 2. ∂H∂uSC=∂H∂uFC in an open interval immediately prior to t∗.
Putting together Lemmas 1 and 2, we get
∂H∂uRC<∂H∂uSC=∂H∂uFC. | (4.15) |
By (3.25), this means that uRC=0, and by (3.38) and (3.39), we obtain
λSG′(zSC)=1 | (4.16) |
during this interval. Next, we state the following result, which means that the shoots are growing during this stage. This will then be used to show that this interval is indeed marked by simultaneous fruit and shoot growth.
Lemma 3. uSN=1 and uSC>0 in an open interval immediately prior to t∗.
At this point, note that we have only showed that
uSC>0,uRC=0,uSN=1,uRN=0. |
Lastly, as λS>1 for this open interval before t∗, λS(t∗)=1, and λS is continuous, we have that
1=limt→t∗−λS=limt→t∗−1G′(zSC)⟹limt→t∗−zSC=0. | (4.17) |
As we have already verified that uSN=1, this means that
limt→t∗−uSC=0, | (4.18) |
and so uSC is continuous at t∗. This also means that 0<uSC<1 during an open interval immediately prior to t∗, so here uFC>0. Therefore, we have established the existence of a penultimate interval, during which
uFC,uSC>0,uRC=0,uSN=1,uRN=0, |
that is, fruits and shoots alone grow simultaneously.
Before discussing the stage that precedes the penultimate interval, we will briefly digress to discuss a key interpretation that will be relevant to the remaining portions of the growth path. One of the most elegant results from the work of Iwasa and Roughgarden in [1] was the relationship between the optimal growth trajectory and the marginal values of the sizes of the organs. In this context, the marginal value of the size of an organ at time t essentially measures the additional fruits present at time T given a unit increase in the size of that organ at time t. In their model, optimal growth was characterized as directing photosynthate to whichever organ had the highest marginal value, and only in the case that two organs had the same marginal value could resources be shared.
Iwasa and Roughgarden identified the adjoints in their model as the marginal values of the shoots and roots, and we can make a similar identification here. Let J(S0,R0,t0) be the maximum for the following optimal control problem.
maxu∫Tt0uFCCdtsubject to:uiC≥0, ujN≥0for i=0,1,2, j=1,2subject to:uFC+uSC+uRC=1=uSN+uRNsubject to:dSdt=g(uSCC,νSuSNN),S(t0)=S0subject to:dRdt=g(uRCC,νRuRNN),R(t0)=R0 | (4.19) |
Note that this only differs from (3.23) in that the initial time is allowed to vary, and so if we initialize (4.19) with a point along the optimal trajectory for (3.23), the optimal trajectories for the two problems will coincide [13,]. In this context, given t∈[0,T] and S(t), R(t) optimal for (3.23), we have that
λS(t)=∂J∂S(S(t),R(t),t) | (4.20) |
λR(t)=∂J∂R(S(t),R(t),t) | (4.21) |
which means that λS(t) and λR(t) measure the increase in the objective function given a small instantaneous increase in either shoots or roots, respectively. See [16] for a more general treatment. Therefore, as J(S(t),R(t),t)=F(T), we can identify λS(t) and λR(t) as the marginal values of the sizes of shoots and roots at time t, respectively. These values change over time, and by (3.31) and (3.37), we know they are both zero at the end of the growing season. This reflects both the compounding effects of investing biomass in shoots/roots over time and the idea that these investments have less of an effect on fruit yield later in the season. Furthermore, as resources invested in fruits are directly converted to tissue, and fruits do not have the capacity to increase their own biomass, we can identify 1 as the marginal value of the size of fruits. That is, a unit of fruits added early on is equivalent to a unit of fruit added near the end of the growing season. It is also worth mentioning that if we changed how we modeled fruits, for example by making fruits photosynthetic, then the marginal value of fruits would change over time to reflect the compounding effects of investing biomass in fruits.
In Iwasa and Roughgarden's model [1], the organ with the highest marginal value was the organ that received the full supply of photosynthate. However, when we consider the optimal trajectory laid out so far, we see that this is not the case in our model, which begs the question why. In particular, in Section 4.2, we see by (4.16) that λSG′(zSC)=1 throughout the interval. Because zSC≠0 on the interior of this interval, we know that G′(zSC)<1 and so λS(t)>1. Furthermore, as λR=0 at the end of the penultimate interval, we have a region where λS(t)>1>λR(t), meaning that the marginal value of the size of shoots is greater than either of the other marginal values, and yet both shoots and fruits are growing.
The reason for this apparent contradiction stems from the addition of the synthesizing units in our model. Consider again the end of the penultimate interval. Here, we see that because the marginal value of shoots is highest, an instantaneous increase in shoots would be better for the plant than either an instantaneous increase in roots or fruits of the same amount. Note, however, that only the control functions can increase instantaneously, and it can be shown that during this phase, the shoot SU alone cannot outperform a mixed shoot/fruit strategy. We will make this more precise below.
Keeping our discussion limited to the penultimate interval, there are two things in particular that we will examine to help us understand why resources are not necessarily directed solely to the organ with the highest marginal value; in this case, shoots. First, we will give a biological interpretation of the necessary conditions during this interval; second, we will use a linear approximation to directly compare the trade-off between a shoot-only strategy and a mixed fruit/shoot strategy during the penultimate interval.
To this end, recall that during the penultimate interval, we have
λSG′(zSC)=1 | (4.16 revisited) |
which relates the marginal values for shoots and fruits. To understand the biological relevance of this statement, we need an interpretation of G′(zSC). By (3.16), we have that
G′(zSC)=∂g∂A(uSCC,νSN)≈g(uSCC+1,νSN)−g(uSCC,νSN) |
during the penultimate interval. That is, G′(zSC) is approximately the additional units of shoot present per unit of time given a unit increase in the rate at which carbon is sent to the shoots. Multiplying this by the marginal value of shoots, λS(t), approximates the additional fruits at time T given the additional shoots produced by a unit increase in carbon allocation to shoots. We can think of the marginal value of fruits being 1 in the same way, in that if we increase the fruit carbon flux, uFCC, by 1, we have an additional
(uFCC+1)−uFCC=1 |
unit of fruit per unit time, which yields an additional unit of fruit at time T. So, (4.16) essentially has the biological meaning that neither the shoot SU nor the fruit SU can outperform the other, and so a mixed strategy is best. In Iwasa and Roughgarden's model [1], raising the rate of carbon flux to an organ by one unit was approximately equivalent in the long run to increasing the biomass of that organ by one unit, which is why comparing marginal values was equivalent to using the necessary conditions.
It is interesting to consider what happens when we attempt to construct a solution based on the marginal values alone during the penultimate interval. To this end, consider a point t# in the penultimate interval such that λS(t#)>1>λR(t#). A growth strategy based solely on the marginal values would suggest shoot-only growth, whereas we know that a mixed fruit/shoot strategy is optimal here based on the necessary conditions. We will examine what happens if at time t# there is a length-ε deviation from the optimal solution during which the mixed fruit/shoot strategy is replaced with a shoot-only strategy. To facilitate this analysis, we have the following theorem.
Theorem 1. Let t# be a point in the penultimate interval such that λS(t#)>1>λR(t#) and such that this ordering of the marginal values holds on [t#,t#+ε] for some ε>0. Let F(T) be the value obtained by following a mixed fruit/shoot strategy on [t#,t#+ε] and then following the optimal trajectory on (t#+ε,T]. Let Falt(T) be the value obtained by following a shoot-only strategy Salt(t) on [t#,t#+ε] and then following the optimal trajectory on (t#+ε,T]. Then, to leading order,
F(T)−Falt(T)≈ε(F′(t#)+λS(t#)(S′(t#)−S′alt(t#)))>0. | (4.22) |
Note that this not only gives us an approximation for the trade-off between the two strategies but also means that to leading order, the solution we have derived is better than the alternative solution suggested solely based on the marginal values. A full proof of Theorem 1 can be found in Appendix B.2, but here we will give an intuitive explanation for the approximation, as well as some numerical evidence.
We can use the marginal values for shoots and fruits at time t# to approximate the impact that each strategy on the interval [t#,t#+ε] will ultimately have on the value of F(T). First, we will consider a mixed fruit/shoot strategy on [t#,t#+ε]; then we have
F(t#+ε)−F(t#)≈εF′(t#)S(t#+ε)−S(t#)≈εS′(t#), |
and so the total contribution to F(T) during the interval [t#,t#+ε] is approximately
εF′(t#)+λS(t#)εS′(t#). |
Likewise, if we consider a shoot-only strategy on [t#,t#+ε], then we have
Salt(t#+ε)−Salt(t#)≈εS′alt(t#), |
and so the total contribution to Falt(T) here is approximately λS(t#)εS′alt(t#). Therefore, assuming that ε is small enough that any difference between the optimal growth paths that the two solutions follow after t=t#+ε is negligible, we conclude that
F(T)−Falt(T)≈ε(F′(t#)+λS(t#)(S′(t#)−S′alt(t#))). |
We use a numerical experiment to illustrate (4.22). Using methods we will discuss in Section 6, we directly compared the difference F(T)−Falt(T) with the approximation ε(F′(t#)+λS(t#)(S′(t#)−S′alt(t#))). Here, we used the terminal conditions S(T)=222.11 and R(T)=111.55 (see Section 7.2). Using ten linearly spaced points in the penultimate interval where λS>1>λR, and starting at each point, we simulated shoot-only growth on a length-ε interval, with ε=10−3, and then continued along the trajectory that is optimal for the new states at the end of each interval. This allowed us to compute F(T)−Falt(T) directly. We also computed the estimate ε(F′(t#)+λS(t#)(S′(t#)−S′alt(t#))) for each of the ten sampled points, and then compared the computed value to the estimate. The results of this experiment are shown in Figure 2. Observe that the estimate follows the computed values, which illustrates the accuracy of our estimate. Note also that F(T)−Falt(T)>0, meaning that the mixed fruit/shoot strategy outperforms the shoot-only strategy, as claimed.
At this point, we have established the existence of a final interval during which only the fruits are growing, and a period before this when both the shoots and fruits are growing together. Continuing backwards we have arrived at a juncture where there are seemingly several options. This is made complicated by the fact that zRC is not defined during the penultimate interval, and so it is not clear directly from the necessary conditions (3.25) which growth pattern is optimal here.
Prior to this period of mixed shoot/fruit growth, there are several possibilities: fruit growth, root/fruit growth, root/shoot/fruit growth, root/shoot growth, shoot-only growth, and root-only growth. Because we assume that fruits are only carbon-dependent, it is never optimal for the plant to invest fully in fruits unless the shoots are done growing. We can also rule out either shoot-only growth or root-only growth here. As we will subsequently see, the first stage of growth consists of either shoot-only growth or root-only growth, and while it may be the case that the length of this second phase is infinitesimally small, we still consider it as distinct from the first phase. That said, it is not optimal to alternate between growing only shoots and only roots, because after addressing any initial imbalance between shoots and roots, to do so would only create further resource limitation for the growing organ.
This leaves us with three remaining options: root/fruit growth, root/shoot/fruit growth, or root/shoot growth. We will use an argument based on our discussion of marginal values in Section 4.3 to eliminate both root/fruit and root/shoot/fruit growth. The necessary conditions (3.25) provide conditions that must be met on an open interval in each case:
root/fruit growth:λRG′(zRC)=1root/shoot/fruit growth:λSG′(zSC)=λRG′(zRC)=1root/shoot growth:λSG′(zSC)=λRG′(zRC)>1. |
As with λSG′(zSC) in Section 4.3, we can think of λRG′(zRC) as approximating the additional units of fruit at time T given the additional roots produced by a unit increase in carbon allocation to roots. When the growth path follows λRG′(zRC)=1, this means that, with respect to carbon, the root SU can do no better for the plant in the long run than the fruit SU can do at present. In the same vein, when a growth path follows λRG′(zRC)>1, the long-term returns of investing carbon in roots exceeds the immediate returns of investing carbon in fruits. This means that root/shoot growth is optimal here, which is consistent with the work of [1].
It is not optimal to have an interval without fruit growth between two intervals that include fruit growth. Such a growth path could be improved by swapping the first two intervals so that the fruits grow continuously. This means that there is no fruit growth prior to the balanced growth phase. As we have already discussed the idea that alternating phases of shoot-only and root-only growth is not optimal, we are left with one option for the initial phase, namely that it consists of either shoot-only growth or root-only growth. Here, the plant addresses deficiencies present in either vegetative organ depending on the initial conditions. Note that this is also the first stage that we see in the single-resource case described in [1].
In this section, we will present the basic equations governing the dynamics in each phase of the solution to (3.23) as well as show that zSC and zRC are continuous between any two consecutive phases in which they are defined. To keep this section organized, we will go through the dynamics for each phase in chronological order and then go on to discuss the transitions between phases. We will also limit this section to include only phase-specific versions of the differential equations for states and adjoints, relevant algebraic constraints due to the necessary conditions, and any additional equations used for understanding the transitions between phases. When we discuss the numerical scheme, we will make use of several additional equations that govern dynamics in various stages but are not relevant here. In what follows, we will also make frequent use of the notation G2(x)=G′(1/x).
During shoot-only growth, we have
uFC=0,uSC=1,uRC=0,uSN=1,uRN=0 | (5.1) |
and only S,λS, and λR are changing, with differential equations in time given by
S′=νSNG(zSC) | (5.2) |
λ′S=−λSCSG′(zSC) | (5.3) |
λ′R=−λSNRνSG2(zSC). | (5.4) |
Here, we have that
zSC=CνSN. | (5.5) |
Note also that in this stage, both N and NR are constant because R is constant. Furthermore, because we know by (4.13) that H=C∗, the Hamiltonian (3.31) during this phase becomes
C∗=λSνSNG(zSC). | (5.6) |
By solving for λS, and employing (3.20), we obtain two expressions for λS during this phase:
λS=C∗νSNG(zSC) | (5.7) |
λS=C∗νSNG2(zSC)+CG′(zSC). | (5.8) |
During root-only growth, we have
uFC=0,uSC=0,uRC=1,uSN=0,uRN=1 | (5.9) |
and only R,λS, and λR are changing, with differential equations in time given by
R′=νRNG(zRC) | (5.10) |
λ′S=−λRCSG′(zRC) | (5.11) |
λ′R=−λRNRνRG2(zRC). | (5.12) |
Here, we have that
zRC=CνRN. | (5.13) |
Note also that in this stage, both C and CS are constant because S is constant. Furthermore, because we know by (4.13) that H=C∗, the Hamiltonian (3.31) during this phase becomes
C∗=λRνRNG(zRC). | (5.14) |
By solving for λR, and employing (3.20), we obtain two expressions for λR during this phase:
λR=C∗νRNG(zRC) | (5.15) |
λR=C∗νRNG2(zRC)+CG′(zRC). | (5.16) |
During balanced growth, we have
uFC=0,0≤uSC≤1,uRC=1−uSC,0≤uSN≤1,uRN=1−uSN | (5.17) |
and so by (3.25), we have that the following two conditions must hold during this phase
λSG′(zSC)=λRG′(zRC)>1 | (5.18) |
λSνSG2(zSC)=λRνRG2(zRC). | (5.19) |
During this stage, S,R,λS, and λR are changing. The differential equations in time for these four during this phase are given by
S′=νSuSNNG(zSC) | (5.20) |
R′=νR(1−uSN)NG(zRC) | (5.21) |
λ′S=−CSλSG′(zSC)=−CSλRG′(zRC) | (5.22) |
λ′R=−NRλSνSG2(zSC)=−NRλRνRG2(zRC). | (5.23) |
Furthermore, again taking advantage of the fact that (4.13) gives us H=C∗, we can rewrite the Hamiltonian (3.31) during this phase as
C∗=λSνSuSNNG(zSC)+λRνRuRNNG(zRC). | (5.24) |
Rewriting this using (3.20), we get
C∗=λSνSuSNN[G2(zSC)+zSCG′(zSC)]+λRνRuRNN[G2(zRC)+zRCG′(zRC)]. |
Simplifying, and using (5.18) and (5.19) to write everything in terms of zSC, gives us
C∗=λSνSuSNN[G2(zSC)+zSCG′(zSC)]+λSuRNN[νSG2(zSC)+νRzRCG′(zSC)]=λS[νSNG2(zSC)(uSN+uRN)+CG′(zSC)(uSC+uRC)]=λS[νSNG2(zSC)+CG′(zSC)]. | (5.25) |
Solving for λS gives us
λS=C∗νSNG2(zSC)+CG′(zSC). | (5.26) |
Using (5.18) and (5.19) to rewrite (5.25) in terms of zRC instead leads to
C∗=λR[νRNG2(zRC)+CG′(zRC)], | (5.27) |
which upon solving for λR gives us
λR=C∗νRNG2(zRC)+CG′(zRC). | (5.28) |
Note that (5.26) and (5.28) are consistent with (5.8) during shoot-only growth, and (5.16) during root-only growth. While we do not assume that either zSC or zRC is continuous between phases (because the controls need not be), we will use this consistency later to show continuity between the initial and balanced growth phases.
During the penultimate interval, we have
0≤uFC≤1,uSC=1−uFC,uRC=0,uSN=1,uRN=0 | (5.29) |
and so by (3.25), we have that
λSG′(zSC)=1. | (4.16 revisited) |
Note that here, we have
zSC=uSCCνSN∗. | (5.30) |
During this interval, S,F,λS, and λR are all changing. The differential equations in time for these four during this stage are
S′=νSN∗G(zSC) | (5.31) |
F′=uFCC | (5.32) |
λ′S=−CS | (5.33) |
λ′R=−N∗RλSνSG2(zSC) | (5.34) |
Furthermore, again taking advantage of the fact that (4.13) gives us H=C∗, we can rewrite the Hamiltonian (3.31) during this phase as
C∗=uFCC+λSνSN∗G(zSC). | (5.35) |
Using (3.20) and (5.30), we get
C∗=uFCC+λSνSN∗(G2(zSC)+zSCG′(zSC))=uFCC+λSνSN∗G2(zSC)+λSG′(zSC)uSCC |
which by (4.16) becomes
C∗=uFCC+uSCC+λSνSN∗G2(zSC)=C+λSνSN∗G2(zSC). | (5.36) |
Now, at this point, we can either solve for λS and obtain
λS=C∗−CνSN∗G2(zSC), | (5.37) |
or we can use (4.16) to rewrite (5.36) as
C∗=CλSG′(zSC)+λSνSN∗G2(zSC) |
and solve to obtain the familiar expression
λS=C∗νSN∗G2(zSC)+CG′(zSC). | (5.38) |
As we commented earlier, the reappearance of this expression for λS suggests that zSC may be continuous between the balanced growth phase and penultimate interval. We will prove this later with the help of (5.37).
Another important result in the penultimate interval, which we prove here, is that zSC is strictly decreasing throughout this phase. This will be especially important because, as we will see in Section 6.3, it will be advantageous for us to think of zSC, rather than t, as the integration variable for the numerical scheme during the penultimate interval. In particular, we will see that numerically solving the differential equations for this phase in t requires integrating a singularity, whereas solving the differential equations in zSC does not.
Lemma 4. zSC is strictly decreasing during the penultimate interval.
Proof. Recall that during the penultimate interval, we have the relationship
λSG′(zSC)=1. | (4.16 revisited) |
Differentiating both sides with respect to t yields
λ′SG′(zSC)+λSG″(zSC)dzSCdt=0. |
Using the fact that λ′S=−CS, we have
dzSCdt=−λ′SG′(zSC)λSG″(zSC)=CS[G′(zSC)]2G″(zSC)=−CS(1+2zSC)26zSC(1+zSC)(1+zSC+z2SC)<0. | (5.39) |
Therefore, zSC is strictly decreasing throughout the penultimate interval.
During the final interval of fruit-only growth, we have
uFC=1,uSC=0,uRC=0,uSN=1,uRN=0 | (5.40) |
and only F and λS are changing, with differential equations in time given by
F′=C∗ | (5.41) |
λ′S=−C∗S, | (5.42) |
Because λS(T)=0, we also have
λS(t)=C∗S⋅(T−t). | (4.6 revisited) |
Lastly, recall that λR(T)=0 and λ′R=0 during this phase, and so λR=0 here as well.
Now that we have established the main dynamics in each phase, we will turn our attention to the transitions between phases. Again, we will proceed in chronological order, beginning with the transition from the initial phase to balanced growth. We can use the fact that the states and adjoints are continuous at this boundary to characterize this transition. In particular, we will show that zSC is continuous at this transition when the first stage consists of shoot-only growth and zRC is continuous at this transition when the first stage consists of root-only growth. In the first case, this will mean that the ratio of the shoot carbon flux to shoot nitrogen flux is continuous across this boundary, and in the second case, the ratio of root carbon flux to root nitrogen flux is continuous.
Due to the symmetry between equations (5.26) and (5.8) and equations (5.28) and (5.16), we will streamline the argument by considering a generalized version of these equations:
λ=C∗νNG2(z)+CG′(z), | (5.43) |
where (λ,ν,z) is either (λS,νS,zSC) or (λR,νR,zRC) depending on whether the initial stage is shoot-only growth or root-only growth, respectively. We call this transition point t=ˉt, let x=limt→ˉt+z, and note that limt→ˉt−z=CνN. Now, taking limits of (5.43) from both sides:
limt→ˉt+λ(t)=ˉλ+=C∗νNG2(x)+CG′(x) | (5.44) |
limt→ˉt−λ(t)=ˉλ−=C∗νNG2(CνN)+CG′(CνN). | (5.45) |
Note that while we do not know the values of C and N at t=ˉt, in what follows it will only matter that they are the same in both (5.44) and (5.45). At t=ˉt, we have ˉλ−=ˉλ+ and so
G2(CνN)+CνNG′(CνN)=G2(x)+CνNG′(x). | (5.46) |
Let
f(x):=G2(x)+CνNG′(x), | (5.47) |
and note that is clear from (5.46) that x=CνN is a solution to
f(x)=G2(CνN)+CνNG′(CνN). | (5.48) |
Here, we will show that this solution is unique. Consider f′(x), recalling (3.21):
f′(x)=G′2(x)+CνNG″(x)=G″(x)(CνN−x). |
Note that because G″(x)≤0 for x≥0, we have that
sgn(f′(x))=sgn(x−CνN)for x≥0. | (5.49) |
Now, as f is continuous and by (5.49), we have that f is decreasing on (0,CνN) and increasing on (CνN,∞), it must be the case that x=CνN is the only solution to (5.48). Therefore, we have that z is continuous at ˉt, which verifies our claim that zSC is continuous at the transition from shoot-only growth to balanced growth and zRC is continuous at the transition from root-only growth to balanced growth.
As in Section 5.6, we can use the fact that the states and adjoints are continuous to show that zSC is continuous at the boundary between the balanced growth and penultimate intervals. Let t=ˆt denote this transition point, and let z−SC and z+SC be the left and right limits of zSC at t=ˆt, respectively. We have, then, the following theorem.
Theorem 2. zSC is continuous at t=ˆt.
We include the details of the proof in Appendix C and only give a sketch of the proof here. To that end, note that by combining (5.26), (5.37), and (4.16), and writing ˆC=C(ˆt), we have that (x,y)=(z−SC,z+SC) must solve
C∗νSN∗G2(y)=(C∗−ˆC)(νSN∗G2(x)+ˆCG′(x)) | (5.50) |
C∗G′(y)=νSN∗G2(x)+ˆCG′(x). | (5.51) |
Because G2 and G′ are both one-to-one, we have that y is a function of x in (5.50) and (5.51).
Figure 3 illustrates the three types of solutions to (5.50) and (5.51), depending on the values of the parameters. Note that the intersection points (x,y) of the two curves correspond to possible pairs (z−SC,z+SC). In particular, there are two common features of each plot:
1. In each plot there is one y value that solves (5.50) and (5.51); this determines z+SC.
2. There are at most two intersection points, one is always at (z+SC,z+SC).
For each case in Figure 3, the only admissible solution is z−SC=z+SC. In Case 2 (Figure 3, center), this is the only possibility. In Case 1 (Figure 3, left), we observe that
ˆCνSN∗<z+SC. |
However, as we know that
z+SC=limt→ˆt+uSCCνSN∗≤ˆCνSN∗, |
this would mean that
ˆCνSN∗<ˆCνSN∗ |
which is impossible. Therefore, Case 1 is not an admissible configuration of solutions to (5.50) and (5.51). In Case 3 (Figure 3, right), there is a possible intersection point at x=p>z+SC. However, if z−SC=p, then as G′ is strictly decreasing, we have that G′(z−SC)<G′(z+SC), and
λS(ˆt)G′(z−SC)<λS(ˆt)G′(z+SC)=1, |
which is impossible by (5.18). Therefore, the only possible solution to (5.50) and (5.51) is z−SC=z+SC, which means that zSC is continuous at t=ˆt. It remains to show that these three cases actually characterize the solution, the details of which are included in Appendix C.
Note that in the final interval, neither zSC nor zRC is defined, so we will not be addressing continuity of either one here. We will, however, analyze the controls at this transition and conclude that at this particular junction, we actually have continuity of the controls.
First, recall that in Section 4.2, we showed that limt→t∗−uSC(t)=0. Because during the penultimate interval uFC=1−uSC we have that limt→t∗−uFC(t)=1 as well. Note also that during the penultimate interval, we have from (5.29) that uRC=0,uSN=1, and uRN=0, and by (5.40) we have that all the controls maintain their values in the final interval as well. Therefore, all the controls are continuous at t=t∗.
Next, we will prove the following lemma, which shows that the transition to fruit-only growth occurs rapidly at the end of the penultimate interval.
Lemma 5.
limt→t∗−d(uSCC)dt=limt→t∗−dzSCdt(t)=−∞. |
Proof. The first equality follows from the penultimate interval definition of zSC (5.30). For the second equality, recall that by (4.17), zSC→0+ as t→t∗−. Lastly, by (5.39), we have
dzSCdt∼−C∗S3z2SC, as zSC→0+, |
which because C∗S>0 means that
limt→t∗−dzSCdt=limzSC→0+−C∗S3z2SC=−∞. |
In this section, we will outline the numerical scheme we developed for solving the optimal control problem (3.23). Because the solution exhibits a multi-phase structure, standard numerical methods such as forward-backward sweep or shooting [13,] are not well-suited. Our approach is to construct a numerical scheme for solving the problem backward in time. There are two primary components to the numerical scheme. First, we obtain the map
(S∗,R∗)↦(S,R,F,u,λS,λR). | (6.1) |
Second, we use MATLAB's built-in nonlinear equation solver fsolve to find the map
(S0,R0)↦(S∗,R∗). |
Ultimately, we obtain the map
(S0,R0)↦(S,R,F,u,λS,λR). |
We will now direct our attention to the numerical scheme for finding (6.1), proceeding phase-by-phase in reverse order. During each phase, we use the fourth-order Runge-Kutta method (RK4) to solve differential equations that govern the phase dynamics and use algebraic equations to update the controls, thus allowing us to avoid having to update the controls iteratively. We will use the fact that the Hamiltonian is constant along the optimal trajectory, as well as the information we have about transitions from Section 5, to find the boundaries between phases and ultimately stitch the four phases together to form one complete solution. We will discuss how the numerical scheme works without getting into the fine numerical details. The MATLAB code is publicly available on David McMorris' GitHub repository, https://github.com/DavidMcMorris/plant-code/blob/main/carbon-only-model.m.
We begin by using equation (4.12) to find the end of the penultimate interval (t∗). Now, as the dynamics during the penultimate interval depend on zSC, we would ideally use (5.39) and (5.31) to solve for zSC and S simultaneously, and use (5.30) to update the controls. However, recall that by Lemma 5, zSC approaches a vertical tangent as t→t∗−, which makes it difficult to solve for zSC numerically. Note, however, that by the same reasoning, we have that dtdzSC approaches 0 as zSC→0+. So, to remedy this numerical difficulty, we take advantage of this fact, along with the fact that, by Lemma 4, zSC is strictly decreasing throughout the penultimate interval, to think of zSC as our ''time" and derive differential equations for t,S, and λR in terms of zSC. This, combined with the fact that during the penultimate interval R is constant, λS is algebraically related to zSC via (4.16), and F depends on the value of ˆt, gives us everything we need to numerically find the solution during this phase. Note that because zSC is decreasing, solving forward in zSC is equivalent to solving backward in time.
To this end, note that by (5.39), we have that
dtdzSC=G″(zSC)CS[G′(zSC)]2,t(0)=t∗. | (6.2) |
Furthermore, using (5.31), (5.34), and (6.2), we obtain
dλRdzSC=dλRdtdtdzSC=−N∗RνSG2(zSC)G″(zSC)CS[G′(zSC)]3,λR(0)=0 | (6.3) |
dSdzSC=dSdtdtdzSC=νSN∗G(zSC)G″(zSC)CS[G′(zSC)]2,S(0)=S∗. | (6.4) |
We use RK4 to solve (6.2), (6.3), and (6.4) forward in zSC, and use (5.30) to update uSC. Also, recall that λS=1/G′(zSC), uFC=1−uSC, and the other three controls are constant by (5.29). Upon reordering by t, this gives us all of the states, adjoints, and controls during the penultimate interval, with the exception of fruits, which we will determine after identifying the correct transition point from balanced growth to the penultimate interval (ˆt).
We stop solving forward in zSC when uSC becomes unbounded, as this gives us an upper bound on the value of zSC (lower bound for the value of t) for which the transition from balanced growth to the penultimate interval can occur. Furthermore, because during the balanced growth phase we have by (5.18) that λRG′(zRC)>1, it must be the case that λR(ˆt)≥1 because G′(z)≤1. Finding where λR=1 gives us a lower bound on the value of zSC (upper bound for the value of t) for which the transition from balanced growth to the penultimate interval can occur. This gives us an interval [zminSC,zmaxSC] in which to search for ˆt.
This portion of the numerical scheme begins with the interval [zminSC,zmaxSC] identified in Section 6.1. Since we know that the Hamiltonian must be constant along the optimal trajectory, we use the balanced-growth specific formulation of H in terms of λR and zRC, given by (5.27):
H=C∗=λR[νRNG2(zRC)+CG′(zRC)]. | (5.27 revisited) |
Starting with [zminSC,zmaxSC], we use a binary search to find the smallest interval, relative to current RK4 step size, which contains the point t=ˆt at which (5.27) is satisfied. Next, we use a smaller step size, and follow the procedure outlined in Section 6.1 to increase the resolution of the solution and repeat the binary search until we find ˆt to within some specified tolerance. This also gives us the values of the states, adjoints, and controls at t=ˆt.
One crucial feature on the process that we have not mentioned so far is the fact that since zRC is not defined in the penultimate interval, we need to compute its limit from balanced growth at every step in the binary search. This is made possible by the fact that zSC is continuous at ˆt, as established in Section 5.7. Therefore, we can use (5.18) and (5.19) to obtain
G′(zRC)G2(zRC)=νRνSG′(zSC)G2(zSC). | (6.5) |
Using MATLAB's built-in nonlinear equation solver fsolve, we use the value of zSC at a particular point to calculate what zRC would be if balanced growth ended at that point. This value of zRC is then used to evaluate the right-hand side of (5.27) in the binary search.
In Section 6.1, we were unable to solve for F because we did not know the value of ˆt. This was resolved in Section 6.2, so at this point, we use RK4 to solve (5.32) forward in time during the penultimate interval. This completes the numerical solution for the penultimate interval.
Before continuing backward in time, we take advantage of the fact that, by the process described in Section 6.3, we now know the value of F at the beginning of the final interval. Here, the controls are constant by (5.40), S,R, and λR are constant, and λS is given by (4.6). Lastly, using the value of F(t∗) obtained as described in Section 6.3, and the fact that F′ is constant during the final interval by (5.41), we have that
F(t)=F(t∗)+C∗⋅(t−t∗). | (6.6) |
Upon locating ˆt as discussed in Section 6.2, we have the value of the states, adjoints, zSC, and zRC at the end of the balanced growth stage. Recall that the controls need not be continuous, so we do not immediately know their values at the end of balanced growth. We can, however, start with the equations for zSC and zRC during balanced growth
zSC=uSCCνSuSNN | (6.7) |
zRC=(1−uSC)CνR(1−uSN)N | (6.8) |
and solve for uSN and uSC:
uSN=CN−νRzRCνSzSC−νRzRC | (6.9) |
uSC=νSzSCNuSNC. | (6.10) |
We derive differential equations for zSC and zRC during this phase, and ultimately use (6.9) and (6.10) to obtain the controls. The differential equations in time for zSC and zRC are given below, and the derivation is included in Appendix D.
dzSCdt=NRνSνRG(zRC)G2(zSC)−CSG′(zSC)[νSG2(zSC)+zRCνRG′(zSC)]G″(zSC)[νSzSC−νRzRC] | (6.11) |
dzRCdt=NRνSG2(zSC)[G′(zRC)]2−CS[G′(zSC)]2G′(zRC)+G″(zSC)G′(zRC)dzSCdtG′(zSC)G″(zRC) | (6.12) |
Beginning at ˆt, we use RK4 to simultaneously solve (5.20), (5.21), (5.22), (5.23), (6.11), and (6.12) backward in time to obtain S,R,λS,λR,zSC, and zRC, respectively. We use (6.10) and (6.9) to eliminate uSC and uSN in (5.20) and (5.21) so that these six differential equations are expressed exclusively in terms of these six variables. We then use (6.10) and (6.9) to obtain uSC and uSN, as well as the rest of the controls via (5.17). We continue until one of the controls exits [0,1], which gives us the earliest time for the beginning of the balanced growth phase.
The next piece of the numerical scheme finds the start of balanced growth, t=ˉt. Recall that in Section 5.6, we showed that either zSC or zRC is continuous at ˉt, depending on which is defined in the initial phase. In either case, the transition must occur at a point where
uSC−uSN=0. | (6.13) |
This is because in the case where the initial phase consists of shoot growth, we have
limt→ˉt+zSC=limt→ˉt+uSCCνSuSNN=limt→ˉt−CνSN⟹limt→ˉt+uSCuSN=1, |
and in the case where the initial phase consists of root growth, we have
limt→ˉt+zRC=limt→ˉt+(1−uSC)CνR(1−uSN)N=limt→ˉt−CνRN⟹limt→ˉt+(1−uSC)(1−uSN)=1. |
We use essentially the same procedure as in Section 6.2 for refining the transition point between the balanced growth and the penultimate intervals, with the exception that finding ˉt does not require computing the limits of any quantities from the earlier phase, as we had to do with zRC at ˆt. We use a binary search to locate the smallest interval about which (6.13) is satisfied, and then use RK4 to solve equations (5.20), (5.21), (5.22), (5.23), (6.11), and (6.12) backward in time on a smaller integration mesh. We repeat this process until we have found a point at which (6.13) is met to within some specified tolerance. We call this point ˉt.
Note that if no such point exists, then the plant begins in balanced growth, in which case there are only three phases. If ˉt exists, the next step in the numerical scheme is to determine the makeup of the initial phase. Using (5.7) and (5.15), we compute the following at ˉt:
|λSνSNG(CνSN)−C∗|and|λRνRNG(CνRN)−C∗|. |
If the first is smaller, then the initial phase consists of shoot growth; if the second is smaller, the initial phase consists of root growth.
If the first phase consists of shoot-only growth, we use RK4 to solve (5.2), (5.3), and (5.4) simultaneously backward in time until we reach t=0. This gives us S,λS, and λR, respectively. The controls are constant here and given by (5.1), and R=R(ˉt).
If the first phase consists of root-only growth, we use RK4 to solve (5.10), (5.11), and (5.12) simultaneously backward in time until we reach t=0. This gives us R,λS, and λR, respectively. The controls are constant here and given by (5.9), and S=S(ˉt).
In this section, we will present some of the primary results from the numerical simulations of the model. In particular, we will look at several ''representative" simulation results that showcase the different strategies a plant can employ to maximize fruit production, as well as help us understand the relationship between initial and terminal conditions. In each simulation, we made the simplifying assumptions that C(S)=S and N(R)=R, which ignore any possibility of self-shading as the plant grows. We chose νR=1 for convenience, and chose νS=1/3, which had the effect of exaggerating the lengths of the various phases, resulting in easier-to-interpret plots. Note that this results in νRνS=3, which is likely a bit higher than data suggests [17,18,]. Lastly, we chose T=10 because in the case that C(S)=S, we have by (4.12) that t∗=9, which again facilitated our interpretation of the numerical results. We will begin with four examples of different strategies that all reach the same optimal value of fruits at time T. The terminal conditions in each example were arbitrarily chosen because they resulted in a solution that illustrated a particular pattern of growth.
Here, we chose terminal conditions S∗=223.20 and R∗=112.86. This results in F(T)=900, S0=17, and R0=84. In this case, we see the full four-stage structure of the solution. There is an initial phase of shoot growth, followed by a period of balanced growth between shoots and roots, a penultimate interval of shoot and fruit growth, and finally a period of fruit-only growth at the end of the growing season. The states and controls are shown in Figure 4.
The lower two plots in Figure 4 show the carbon and nitrogen allocation strategies for this plant. Following the initial period of shoot growth, there is a short period during which uRC and uRN are increasing, signifying a period of increasing root production. Shortly after t=2, we see uRC and uRN begin to decrease, signifying that although it is still advantageous to invest resources in roots, the plant ultimately needs to prioritize shoot production again to prepare for mixed shoot/fruit growth. During the penultimate interval, we see the plant gradually stop investing in shoots, before switching to fruit-only growth at time t=9.
Here, we chose terminal conditions S∗=222.11 and R∗=111.55. This results in F(T)=900, S0=48.9, and R0=28.9. In this simulation, we again see all four stages; however, here we see an initial phase of root growth instead of the initial phase of shoot growth we saw previously. The initial conditions of the two simulations are quite different, but the terminal conditions are nearly the same. This, as we will discuss later, suggests that optimal growth tends to smooth out initial transients. The states and controls are shown in Figure 5. The two graphs at the bottom of Figure 5 show a different allocation strategy than what appeared in the case of initial shoot growth in Figure 4. In particular, following the initial period of root growth, we see a decline in both carbon and nitrogen allocation to the roots throughout the entire balanced growth phase, and the steady increase in allocation to shoots.
Here, we chose terminal conditions S∗=222.71 and R∗=112.27. This results in F(T)=900, S0=30.3, and R0=59.74. The states and controls are shown in Figure 6. In this simulation, we see a plant that starts in balance and therefore forgoes the initial phase of shoot-only or root-only growth. That said, comparing Figures 6 and 4, we see a similar allocation strategy during balanced growth. In particular, both show that allocation to shoots initially decreases before increasing again. This behavior corresponds to initial conditions that are more biased toward shoot deficiency than root deficiency. For this reason, we refer to this type of initially balanced growth as Type S.
Here, we chose terminal conditions S∗=222.60 and R∗=112.14. This results in F(T)=900, S0=36.16, and R0=49.58. The states and controls are shown in Figure 7. Here, we see another example of a plant that begins in balanced growth. Unlike Figure 6, the balanced growth phase seen here has a similar structure to that in Figure 5, where the plant begins with root-only growth. In both examples, the phase consists of steadily increasing allocation to shoots while simultaneously decreasing allocation to roots. In this case, the plant is initially more biased toward root deficiency than shoot deficiency and so the early growth sees a greater investment in roots. We refer to this type of initially balanced growth as Type R.
In order to better understand the relationship between initial conditions and optimal fruit yield, we looked for points in the (S0,R0)-plane for which F(T) is 700,800, or 900. For a given value of S0, we used MATLAB to find the corresponding value of R0 for which the numerical scheme outlined in Section 6 would yield either 700,800, or 900 for F(T). We plotted the resulting contours in the (S0,R0)-plane as seen in Figure 8. We also note that the four representative simulations we have examined occur on the F(T)=900 contour.
There are several key features of this plot to notice. First is the wide range of initial conditions for which a particular fruit yield is optimal. Depending on the allocation strategy, plants with very different initial conditions may share an optimal fruit yield. Additionally, for any given point in Figure 8, we can see how much more initial shoots or roots would be necessary to move the plant to a higher-yield contour. Looking at the ends of the contours, it takes relatively little additional initial structure to move from one contour to the next, whereas in the middle of the contours we see that a larger increase in initial biomass is required to move to the next contour. Looking at Figures 6 and 7, we see that these middle regions contain initial conditions from which the plant begins in balance. Nearer to the ends of the contours, where the concavity becomes more pronounced, we see that the plant must be more biased toward the vegetative organ with the most biomass. It is also worth noting that due to the level of generality with which we built the model, it is difficult to say where the boundaries lie in the (S0,R0) plane for realistic initial conditions, and it may be biologically unlikely for a plant to begin in balance.
By incorporating nitrogen into the model, but keeping the fruits solely reliant on carbon, we have considered a case in which νF→∞. While not realistic for most annual plants, this provides a mathematically approachable framework in which to begin to analyze optimal allocation of two resources. It is worth noting, however, that if we exchange the assumption that fruits are C-only with the assumption that the fruits are N-only, we expect the model to predict a reversal in the roles of shoots and roots. The assumption that fruits are C-only, however, is a more natural extension of Iwasa and Roughgarden's work [1].
An obvious outcome of our model, not present in [1], is the penultimate interval of mixed shoot/fruit growth. As νF>νS,νR, we see a phase during which shoot production is more important overall to eventual fruit yield, but an increased capacity to assimilate nitrogen is not useful, so any excess carbon is invested in fruits. As we discussed in Section 4.3, this corresponds to when the marginal value of the size of shoots exceeds the marginal values of the sizes of the other organs, but because the plant allocates resources to synthesizing units in each organ rather than biomass directly, we see a situation where allocating carbon to both shoots and fruits is more beneficial than allocating only to shoots. We also note that this behavior is consistent with recent predictions in the literature [19].
The theoretical results in Section 5 confer a degree of biological relevance to the model that adds credence to the model despite the narrow scope imposed by the assumption that fruits are built solely from carbon. In particular, we showed that zSC,zRC, and uFC are continuous between any two phases in which they are defined and nonzero. Note that while zSC and zRC are dimensionless, they are multiples of the ratio of C flux to N flux in shoots and roots. So, the fact that zSC and zRC are continuous here means that these ratios are continuous in the shoots and roots. This is striking, given that, while these individual fluxes are continuous at the beginning of balanced growth, they are markedly discontinuous at the end in Figures 4, 5, 6, and 7. So, as long as either shoots or roots are growing, the amount of allocated C per unit of allocated N varies continuously, which is reasonable to expect from biochemical processes. Furthermore, by (3.30), the fact that uFC is continuous between the penultimate and final intervals means that the rate of fruit growth is continuous here.
The results presented in Section 7 provide several avenues for drawing more general conclusions about the nature of plant growth that optimizes fruit yield. In Figures 4, 5, 6, and 7, we see that allocation is a balancing act between avoiding limitation and investing in the organ that most directly contributes to increased fruit yield. The plant invests in the most deficient organ first until it can efficiently invest in both shoots and roots, doing so that by the end of balanced growth, we see increased C flux to shoots. During the penultimate interval, fruits are a better investment than roots, but fruits still benefit from increased C flux. Therefore, the plant continues to invest in shoots while gradually transitioning to fruit-only growth. These ideas are also reminiscent of our discussion of the marginal values and necessary conditions for optimality in Section 4.3, in that they suggest that two organs can grow simultaneously only if neither SU can outperform the other.
Recall that Figures 4, 5, 6, and 7 represent snapshots of optimal growth along the F(T)=900 contour represented in Figure 8. In particular, we carefully chose initial conditions and parameter values to present these snapshots with clearly-defined phases. However, it is worth noting that in general, the length of each phase can vary. It is also plausible that given initial conditions that are sufficiently biased toward either shoots or roots, along with a short enough growing season, the plant may jump from the initial stage directly to the final phase of reproductive growth, in which case we would have bang-bang controls. This, however, is not something we have analyzed from a mathematical standpoint, as in practice, a phase being omitted is biologically equivalent to it being infinitesimally short due to the timescales involved with resource allocation. It is also unclear how successful a plant could be in such an environment where the length of the growing season precludes the possibility of balanced growth. That is to say, while the general trends discussed above hold across the contour, there is a broad spectrum of optimal strategies. What is striking is the high level of variation in solutions with initial conditions along the contour, and the fact that these initial transients are absent by time T. We see this in the values of S(T), R(T), and F(T) in Figures 4, 5, 6, and 7, but the trend generally holds along the contour. In some sense, then, we can think of optimal growth under this model as being an equalizing agent that reduces initial variation in a population.
A natural extension to this model is to drop the assumption that fruits are carbon-only and consider the more general case where fruits require both carbon and nitrogen. Here, we expect the optimal growth path to depend on the stoichiometric ratios of carbon to nitrogen, with a penultimate interval consisting of either shoot/fruit or root/fruit growth depending on the relative importance of each resource to the fruits.
The inclusion of nitrogen also opens other avenues for inquiry that have been previously explored in carbon-only optimal allocation models. For example, we have not included the possibility of resource storage and remobilization in our model. Chiariello and Roughgarden studied a carbon-only allocation model for annual plants in [20] and found that storing carbon for future allocation to reproduction can be optimal in environments where the rate of seed loss exceeds certain thresholds related to the net growth rate of the structural tissue, or when the efficiency with which carbon can be stored/remobilized is sufficiently high. It seems likely that a similar conclusion would hold with the inclusion of nitrogen, although the conditions under which storage of carbon, nitrogen, or both is optimal would likely depend on the resource requirements for building fruits as well as how efficiently the plant can store/remobilize each resource. Investment in storage may also be optimal in an environment in which the plant may experience sudden tissue loss due to disturbances, in which case the storage compartment acts as a safety net to allow the plant to recover [21]. Here, again, we would generally expect the inclusion of a second resource in the model to allow for richer storage dynamics, where the size of the reserves may depend on the resource requirements for rebuilding lost tissue and the frequency/severity of disturbances.
In this appendix, we will derive the necessary conditions presented in (3.25) in Section 3.4. Note that because the carbon controls do not directly depend on the nitrogen controls, we can obtain the necessary conditions by considering an n-state two-control problem and an n-state three-control problem separately. In our particular situation, n=3, but we derive the conditions with higher generality because there is minimal complexity added to the derivation. In the following, we consider f,g to be continuously differentiable in all arguments and u to be piecewise continuous. We begin with an n-state two-control problem in Appendix A.1 followed by an n-state three-control problem in Appendix A.2.
We consider the optimal control problem given by (A.1), with the goal of deriving the associated necessary conditions.
maxu=(u1,u2)∫T0f(t,x(t),u(t))dtsubject to:0≤u1(t)≤1,0≤u2(t)≤1for t∈[0,T],subject to:u1(t)+u2(t)=1for t∈[0,T]subject to:x′(t)=g(t,x(t),u(t)),x(0)=x0 | (A.1) |
To this end, let U be the space of admissible controls, and define the functional J:U→R by
J(u)=∫T0f(t,x(t),u(t))dt=∫T0{f(t,x(t),u(t))+λ(t)⋅(g(t,x(t),u(t))−x′(t))}dt=∫T0{f(t,x(t),u(t))+λ(t)⋅g(t,x(t),u(t))+λ′(t)⋅x(t)}dt=∫T0{f(t,x(t),u(t))+λ(0)⋅x0−λ(T)⋅x(t), | (A.2) |
where λ(t) is a piecewise differentiable function to be specified later, and the final equality was obtained by integrating by parts. Note that by λ′ we mean the function whose nth component is the time derivative of the nth component of λ. We use the same notation for x′(t). Let u∗,x∗ be optimal, and for piecewise continuous variations h1,h2 and ε∈R define (uε1,uε2)∈U by uε1(t)=u∗1(t)+εh1(t) and uε2(t)=u∗2(t)+εh2(t), and let xε be the corresponding state. Note that admissibility requires uε1+uε2=1, so we have h1=−h2. Because u∗,x∗ are optimal, we have
0≥dJ(uε)dε|ε=0=limε→0J(uε)−J(u∗)ε. | (A.3) |
We have an inequality rather than equality above because the boundedness of the controls makes it possible that optimality is only global rather than local, and the numerator in the difference quotient is nonpositive because J(u∗) is maximal over admissible controls. We differentiate (A.2) with respect to ε and note that conditions on the functions involved in the integrand allow us to interchange the derivative and integral via the dominated convergence theorem. Suppressing the arguments from here on,
0≥∫T0∂∂ε(f+λ⋅g+λ′⋅xε)|ε=0dt−λ(T)⋅∂xε(t)∂ε|ε=0=∫T0{(∇f+n∑i=1λi∇gi+λ′)⋅∂xε(t)∂ε|ε=0=∫T0{(∇f++(fu1+λ⋅gu1)h1+(fu2+λ⋅gu2)h2n∑i=1}dt=∫T0{(∇f+(fu1+−λ(T)⋅∂xε(t)∂ε|ε=0. |
Note that here we are using ∇f to mean the vector of derivatives of f with respect to each component of x, and likewise for ∇gi. Choosing λ to satisfy
λ′(t)=−(∇f+n∑i=1λi∇gi),λ(T)=0, |
the inequality reduces to
0≥∫T0{(fu1+λ⋅gu1)h1+(fu2+λ⋅gu2)h2}dt. |
Using the fact that h1=−h2, we obtain
0≥∫T0(fu1+λ⋅gu1−fu2−λ⋅gu2)h1dt. |
Next, we use this inequality to obtain necessary conditions for optimality. Let s be a point of continuity for u∗1 and u∗2 such that 0≤u∗1(s)<1, and suppose for the sake of contradiction that fu1+λ⋅gu1−fu2−λ⋅gu2>0 at s. As λ is continuous, and f and g are continuously differentiable, by the continuity of u∗1 and u∗2 at s, we have that fu1+λ⋅gu1−fu2−λ⋅gu2 is continuous at s as well. Therefore, we can find a closed interval I about s such that u∗1<1 and fu1+λ⋅gu1−fu2−λ⋅gu2>0 on I. Let
M=max{u∗1(t)∣t∈I}<1, |
and choose h1 and h2 to be
h1(t)=(1−M)χI(t),h2(t)=(M−1)χI(t), |
where χI is the characteristic function on I. Note that this gives us
uε1=u∗1+ε(1−M)χI,uε2=u∗2+ε(M−1)χI. |
As u∗1≤M on I, then it must be the case that u∗2≥1−M on I since u∗2=1−u∗1. This means that for all ε∈[0,1] we have
0≤uε1≤1,0≤uε2≤1, |
and as the variations were chosen such that h1+h2=0, we have that these variations lead to admissible controls. We have then
0≥∫T0(fu1+λ⋅gu1−fu2−λ⋅gu2)h1dt=∫I(fu1+λ⋅gu1−fu2−λ⋅gu2)(1−M)dt>0, |
a contradiction. So, it must be the case that fu1+λ⋅gu1−fu2−λ⋅gu2≤0. Note that because the controls sum to one, we have that 0<u∗2≤1 implies 0≤u∗1<1, and so we actually have the following condition
0≤u∗1<1or0<u∗2≤1⟹fu1+λ⋅gu1−fu2−λ⋅gu2≤0. |
Interchanging the roles of u1 and u2 in the argument above also gives us
0<u∗1≤1or0≤u∗2<1⟹fu1+λ⋅gu1−fu2−λ⋅gu2≥0. |
So, together we have that
{0≤u∗1<1or0<u∗2≤1⟹fu1+λ⋅gu1−fu2−λ⋅gu2≤00<u∗1≤1or0≤u∗2<1⟹fu1+λ⋅gu1−fu2−λ⋅gu2≥0 | (A.4) |
Now, suppose that we are in the case where fu1+λ⋅gu1−fu2−λ⋅gu2<0. By (A.4) we can rule out the possibility that either u∗1>0 or u∗2<1 because either choice would lead to a contradiction. Therefore, in this case, we can conclude that u∗1=0 and u∗2=1. Similar arguments lead to the following conditions:
{fu1+λ⋅gu1−fu2−λ⋅gu2<0⟹u∗1=0,u∗2=1fu1+λ⋅gu1−fu2−λ⋅gu2>0⟹u∗1=1,u∗2=0fu1+λ⋅gu1−fu2−λ⋅gu2=0⟹0<u∗1,u∗2<1. | (A.5) |
We can convert this into conditions involving a Hamiltonian as follows. Define H by
H(t,x,u,λ):=f(t,x,u)+λ⋅g(t,x,u). |
For distinct indices i,j∈{1,2}, rewriting (A.5) in terms of H yields the necessary conditions
{∂H∂ui<∂H∂uj⟹u∗i=0,u∗j=1∂H∂ui=∂H∂uj⟹0≤u∗i,u∗j≤1. | (A.6) |
We can also express the differential equations for x and λ in terms of the Hamiltonian as
{x′i(t)=∂H∂λi,xi(0)=xi0 for i=1,…,nλ′i(t)=−∂H∂xi, λi(T)=0 for i=1,…,n. | (A.7) |
Note that if we were to reduce the problem to one control by writing u2=1−u1, these conditions would become exactly the standard necessary conditions for a problem with one bounded control and n states [13,]. Going through this derivation with two controls, however, gives us a starting point to approach the analogous problem with three controls.
We now consider the following optimal control problem with n states and three controls. We will again derive the necessary conditions, following a similar procedure to that used in the simpler case with only two controls. The problem is stated as follows.
maxu=(u1,u2,u3)∫T0f(t,x(t),u(t))dtsubject to:0≤u1(t)≤1,0≤u2(t)≤1,0≤u3(t)≤1for t∈[0,T],subject to:u1(t)+u2(t)+u3(t)=1for t∈[0,T]subject to:x′(t)=g(t,x(t),u(t)),x(0)=x0 | (A.8) |
The derivation of the necessary conditions is similar to the two-control case. We let U be the space of admissible controls for (A.8) and define J:U→R by (A.2), restated below.
J(u)=∫T0{f(t,x(t),u(t))+λ(t)⋅g(t,x(t),u(t))+λ′(t)⋅x(t)}dt=∫T0{f(t,x(t),u(t))+λ(0)⋅x0−λ(T)⋅x(t). | (A.2 revisited) |
Suppose that u∗,x∗ are optimal, and let h1,h2, and h3 be piecewise continuous variations. Then, for ε∈R, we define uε∈U by uεi(t)=u∗i(t)+εhi(t) for i=1,2,3, and let xε be the corresponding state. Admissibility here requires that ∑3i=1uεi=1, so we have that h3=−(h1+h2). By the optimality of u∗,x∗ we again get (A.3), restated below:
0≥dJ(uε)dε|ε=0=limε→0J(uε)−J(u∗)ε. | (A.3 revisited) |
We again differentiate (A.2), using the DCT to interchange the order of differentiation and integration, this time arriving at
0≥∫T0∂∂ε(f+λ⋅g+λ′⋅xε)|ε=0dt−λ(T)⋅∂xε(t)∂ε|ε=0=∫T0{(∇f+n∑i=1λi∇gi+λ′)⋅∂xε(t)∂ε|ε=0+(fu1+λ⋅gu1)h1=∫T0{(∇f++(fu2+λ⋅gu2)h2+(fu3+λ⋅gu3)h3n∑i=1}dt−λ(T)⋅∂xε(t)∂ε|ε=0. |
As before, we choose λ so that
λ′(t)=−(∇f+n∑i=1λi∇gi),λ(T)=0, |
and so we arrive at the inequality
0≥∫T0{(fu1+λ⋅gu1)h1+(fu2+λ⋅gu2)h2+(fu3+λ⋅gu3)h3}dt. | (A.9) |
We will use this inequality to derive the necessary conditions for optimality. Due to the similarity between the various cases, we will only show one in detail. We begin by using the substitution h3=−(h1+h2) to rewrite (A.9) as
0≥∫T0{(fu1+λ⋅gu1−fu3−λ⋅gu3)h1+(fu2+λ⋅gu2−fu3−λ⋅gu3)h2}dt. | (A.10) |
Next, let s be a point of continuity for all controls so that 0≤u∗1(s)<1 and 0<u∗3(s)≤1. Note that because the controls sum to one, this also means that 0≤u∗2(s)<1. Additionally, assume for the sake of contradiction that fu1+λ⋅gu1−fu3−λ⋅gu3>0 at s. As λ is continuous, and f and g are continuously differentiable, by the continuity of the controls, we have that fu1+λ⋅gu1−fu3−λ⋅gu3 is continuous at s as well. Therefore, we can find a closed interval I about s such that u∗1<1,u∗3>0, and fu1+λ⋅gu1−fu3−λ⋅gu3>0 on I. Next, let
m=min{1−max{u∗1(t)∣t∈I},min{u∗3(t)∣t∈I}}. |
and note that m>0. Next, we choose variations to be
h1(t)=mχI(t),h2(t)≡0,h3(t)=−mχI(t). |
This gives us the controls
uε1(t)=u∗1(t)+εmχI(t),uε2(t)=u∗2(t),uε3(t)=u∗3(t)−εmχI(t). |
Note that because ∑3i=1u∗i(t)=1 and ∑3i=1hi(t)=0, we have that ∑3iuεi(t)=1 as well. Furthermore, restricting our attention to t∈I and ε∈[0,1], we have
0≤uε1(t)=u∗1(t)+εm≤u∗1(t)+ε(1−max{u∗1(t)})≤u∗1(t)+1−max{u∗1(t)}≤1. |
As this bound also holds outside of I by assumption, we have that 0≤uε1≤1. Now, for uε3, again restricting our attention to t∈I and ε∈[0,1], we have
1≥uε3(t)=u∗3(t)−εmχI(t)≥u∗3(t)−εmin{u∗3(t)}≥u∗3(t)−min{u∗3(t)}≥0. |
Again, as this bound holds outside of I by assumption, we have that 0≤uε3≤1. Therefore, the controls uε1,uε2, and uε3 are admissible for all ε∈[0,1]. We have then
0≥∫T0{(fu1+λ⋅gu1−fu3−λ⋅gu3)h1+(fu2+λ⋅gu2−fu3−λ⋅gu3)h2}dt=∫I(fu1+λ⋅gu1−fu3−λ⋅gu3)mdt>0, |
a contradiction. Therefore, it must be the case that fu1+λ⋅gu1−fu3−λ⋅gu3≤0.
Now, by swapping the indices 1 and 2 in the argument above, we can reach the additional conclusion that fu2+λ⋅gu2−fu3−λ⋅gu3≤0 for the same set of assumptions. Note that this is because if we assume that 0<u∗3≤1, then we can conclude both 0≤u∗1<1 and 0≤u∗2<1 regardless of which assumption is used for a particular argument. Therefore, if we are in the case that 0≤u∗1<1, 0≤u∗2<1, and 0<u∗3≤1, then we have that both fu1+λ⋅gu1−fu3−λ⋅gu3≤0 and fu2+λ⋅gu2−fu3−λ⋅gu3≤0. We can permute the roles of u1,u2, and u3, or, equivalently, repeat the argument above with the substitutions of h1=−(h3+h2) or h2=−(h3+h1) into (A.9) to obtain the following.
{0≤u∗1,u∗2<1,0<u∗3≤1⟹{fu1+λ⋅gu1−fu3−λ⋅gu3≤0fu2+λ⋅gu2−fu3−λ⋅gu3≤00≤u∗1,u∗3<1,0<u∗2≤1⟹{fu1+λ⋅gu1−fu2−λ⋅gu2≤0fu3+λ⋅gu3−fu2−λ⋅gu2≤00≤u∗2,u∗3<1,0<u∗1≤1⟹{fu2+λ⋅gu2−fu1−λ⋅gu1≤0fu3+λ⋅gu3−fu1−λ⋅gu1≤0 | (A.11) |
We will now work to develop implications going in the other direction. We will illustrate this for a few cases, and the remaining cases follow from permuting the controls. First, consider the case that at a point t∈[0,T], we have fu1+λ⋅gu1−fu3−λ⋅gu3<0 and fu2+λ⋅gu2−fu3−λ⋅gu3<0, and suppose for the sake of contradiction that u∗3(t)<1. Then it must be the case that at t either 0<u∗1≤1 and 0≤u∗2<1 or 0≤u∗1<1 and 0<u∗2≤1. Using (A.11), we see that the first case implies that fu3+λ⋅gu3−fu1−λ⋅gu1≤0 at t, and the second case implies that fu3+λ⋅gu3−fu2−λ⋅gu2≤0 at t, both of which contradict our assumptions. So, it must be the case that u∗3(t)=1.
Next, consider the case that at a point t∈[0,T], we have fu1+λ⋅gu1−fu3−λ⋅gu3<0 and fu1+λ⋅gu1−fu2−λ⋅gu2<0, and suppose for the sake of contradiction that u∗1(t)>0. Then it must be the case that at t, we have 0≤u∗2<1 and 0≤u∗3<1. Using (A.11), we see that this implies that fu2+λ⋅gu2−fu1−λ⋅gu1≤0 and fu3+λ⋅gu3−fu1−λ⋅gu1≤0 at t, both of which contradict our assumptions here. Therefore, it must be the case that u∗1(t)=0. A similar argument can be used to show that if fu2+λ⋅gu2−fu1−λ⋅gu1<0 and fu2+λ⋅gu2−fu3−λ⋅gu3<0 at a point t∈[0,T], then it must be the case that u∗2(t)=0.
By permuting the roles of the controls in the arguments above, we can reach the following conclusions. For distinct indices i,j,k∈{1,2,3}, we have that
{{fuj+λ⋅guj−fui−λ⋅gui<0fuk+λ⋅guk−fui−λ⋅gui<0⟹u∗i=1{fui+λ⋅gui−fuj−λ⋅guj<0fui+λ⋅gui−fuk−λ⋅guk<0⟹u∗i=0 | (A.12) |
Note that this also tells us that the only way for a control to take on a value other than zero or one is for at least one of the inequalities in (A.12) to be an equality instead. We can convert this into conditions involving a Hamiltonian as follows. Define H by
H(t,x,u,λ):=f(t,x,u)+λ⋅g(t,x,u). |
Rewriting (A.12) in terms of H yields the necessary conditions
{∂H∂uj<∂H∂uiand∂H∂uk<∂H∂ui⟹u∗i=1∂H∂ui<∂H∂ujand∂H∂ui<∂H∂uk⟹u∗i=0∂H∂ui=∂H∂ujor∂H∂ui=∂H∂ukor∂H∂uj=∂H∂uk⟹0≤u∗i≤1. | (A.13) |
We can also express the differential equations for x and λ in terms of the Hamiltonian:
{x′i(t)=∂H∂λi,xi(0)=xi0 for i=1,…,nλ′i(t)=−∂H∂xi, λi(T)=0 for i=1,…,n. | (A.14) |
In this appendix, we include the proofs of Lemma 1, Lemma 2, Lemma 3, and Theorem 1.
Lemma 1. ∂H∂uRC<∂H∂uFC in an open interval immediately prior to t∗.
Proof. Because λR(t∗)=0, λR is continuous, and G′ is bounded, there exists ε>0, such that for all t in (t∗−ε,t∗), we have
λRG′(zRC)C<C, |
which in terms of the Hamiltonian (3.31) gives the desired result:
∂H∂uRC<∂H∂uFC. |
Lemma 2. ∂H∂uSC=∂H∂uFC in an open interval immediately prior to t∗.
Proof. We will begin by showing that
∂H∂uSC≥∂H∂uFC |
in an open interval immediately prior to t∗. Suppose to the contrary that for some ε>0 we have for all t in (t∗−ε,t∗) that
∂H∂uSC<∂H∂uFC. |
By Lemma 1, we have that
∂H∂uRC<∂H∂uFC |
as well, and so by (3.25), this means that uFC=1 and uSC=0=uRC. As in the final interval, this yields
∂H∂uFC=C,∂H∂uSC=λSC. |
However, as λS(t∗)=1 and by (3.45), we have that λ′S<0, it must be the case that λS>1. This then means that
∂H∂uSC>∂H∂uFC, |
a contradiction. Therefore, we have shown that during this interval, we have
∂H∂uSC≥∂H∂uFC. |
Suppose now that the inequality is strict, and that in fact we have
∂H∂uSC>∂H∂uFC. |
Now, combining this supposition with Lemma 1 gives us
∂H∂uSC>∂H∂uFC>∂H∂uRC, |
which by (3.25) means that
uSC=1,uFC=0=uRC. |
This then gives us
zSN=νSuSNNC, |
which, for fixed C and N, is bounded between 0 and νSNC. Therefore, we have that G′(zSN)≠0. Furthermore, recall that we know that λR(t∗)=0 and G′ is bounded. This means that there is a (potentially smaller) open interval immediately prior to t∗ in which we have
λSνSNG′(zSN)>λRνRNG′(zRN), |
or in terms of the Hamiltonian,
∂H∂uSN>∂H∂uRN. |
By (3.25), this means that uSN=1 and uRN=0 during this interval. With this configuration of the controls, we have that N=N∗ is constant because R′=0 and this interval reaches this final interval. So, we have that
zSC=CνSN∗>0. |
This then means that G′(zSC)<1, and so as λS(t∗)=1 and both zSC and λS are continuous here, there exists ε>0 such that for all t in (t∗−ε,t∗), we have
λSG′(zSC)<1 |
that is
∂H∂uSC<∂H∂uFC, |
a contradiction. Therefore, we have
∂H∂uSC=∂H∂uFC |
in an open interval immediately prior to t∗, as desired.
Lemma 3. uSN=1 and uSC>0 in an open interval immediately prior to t∗.
Proof. We will first show that uSN is nonzero. Recall that by (4.15) and (3.25), we have that uRC=0 in an open interval immediately prior to t∗. Suppose furthermore that uSN=0 during this interval as well. By (3.6) and (3.7), this means that both S′=0 and R′=0. As this interval is followed by a period of fruit-only growth, in order to maximize the integral in (3.23), it must be the case that uFC=1. We therefore have
uFC=1,uSC=0=uRC,uSN=0,uRN=1. |
In particular, because here we are again in a case where uFC=1, uSC=0=uRC, and λR(t∗)=0, we again arrive at a situation where
∂H∂uFC=C,∂H∂uSC=λSC |
following the same argument as in the justification for the existence of a final interval of fruit-only growth. However, as we know that here λS>1, we have that
∂H∂uFC<∂H∂uSC, |
which, regardless of the relationship between ∂H∂uRC with the other two partial derivatives, contradicts the fact that uFC=1. Hence, uSN>0.
Now, because λS>1, we know by (4.16) that zSC>0 so that G′(zSC)<1. Because uSN>0, this allows us to conclude that uSC>0 as well without any concern regarding ambiguous limiting behavior with G′(zSC) when one or both controls are zero. This in turn implies that
∂H∂uSN=λSνSNG′(zSN)>0. |
Now, if uSN<1, we would also have uRN>0. Because uRC=0, this would imply that zRN→∞, and so by (3.18), we would have G′(zRN)=0. This would then imply that
∂H∂uSN>∂H∂uRN, |
which by (3.25) means that uSN=1, a contradiction. Therefore, it must have been the case that uSN=1 in addition to uSC>0 in an open interval immediately prior to t∗, as desired.
Theorem 1. Let t# be a point in the penultimate interval such that λS(t#)>1>λR(t#) and such that this ordering of the marginal values holds on [t#,t#+ε] for some ε>0. Let F(T) be the value obtained by following a mixed fruit/shoot strategy on [t#,t#+ε] and then following the optimal trajectory on (t#+ε,T]. Let Falt(T) be the value obtained by following a shoot-only strategy Salt(t) on [t#,t#+ε] and then following the optimal trajectory on (t#+ε,T]. Then, to leading order,
F(T)−Falt(T)≈ε(F′(t#)+λS(t#)(S′(t#)−S′alt(t#)))>0. | (4.22 revisited) |
Proof. To derive the approximation in Theorem 1, we will first recall a definition we introduced in Section 4.3 when discussing the marginal values. Let J(S0,R0,t0) be the maximum for the following optimal control problem.
maxu∫Tt0uFCCdtsubject to:uiC≥0, ujN≥0 for i=0,1,2, j=1,2subject to:uFC+uSC+uRC=1=uSN+uRNsubject to:dSdt=g(uSCC,νSuSNN),S(t0)=S0subject to:dRdt=g(uRCC,νRuRNN),R(t0)=R0 | (4.19 revisited) |
We wish to compute F(T)−Falt(T), which in terms of J is given by
∫t#+εt#F′(t)dt+J(S(t#+ε),R∗,t#+ε)−J(Salt(t#+ε),R∗,t#+ε), | (B.1) |
where we use the 'alt' subscript wherever we need to distinguish variables that differ from the optimal solution based on the shoot-only growth strategy during [t#,t#+ε]. Note also that here R=R∗, because we are in the penultimate interval. We can interpret (B.1) as the sum of the amount of fruit grown during the interval [t#,t#+ε] using the mixed fruit/shoot strategy and the difference between the final value of fruits depending on the different values of S(t#) produced by each strategy on [t#,t#+ε]. Now, for small enough ε, we have that
S(t#+ε)≈S(t#)+εS′(t#)=S(t#)+ενSR∗G(uSC(t#)C(S(t#))νSR∗)Salt(t#+ε)≈S(t#)+εS′alt(t#)=S(t#)+ενSR∗G(C(S(t#))νSR∗) |
which makes (B.1) approximately
∫t#+εt#F′(t)dt+J(S(t#)+ενSR∗G(uSC(t#)C(S(t#))νSR∗),R∗,t#+ε)−J(S(t#)+ενSR∗G(C(S(t#))νSR∗),R∗,t#+ε). | (B.2) |
We will simplify the second two terms in (B.2), using the following labels for convenience:
I1=∫t#+εt#F′(t)dtI2=J(S(t#)+ενSR∗G(uSC(t#)C(S(t#))νSR∗),R∗,t#+ε)I3=J(S(t#)+ενSR∗G(C(S(t#))νSR∗),R∗,t#+ε). |
Beginning with I2, and integrating by parts:
I2=∫Tt#+εuFC(t)C(S(t))dt=∫Tt#+ε{uFC(t)C(S(t))+λS(t)S′(t)−λS(t)S′(t)+λR(t)R′(t)−λR(t)R′(t)}dt=∫Tt#+ε{uFC(t)C(S(t))+λS(t)S′(t)+λ′S(t)S(t)+λR(t)R′(t)+λ′R(t)R(t)}dt−λS(T)S(T)+λS(t#+ε)S(t#+ε)−λR(T)R(T)+λR(t#+ε)R(t#+ε). |
Likewise for I3, using the same λS,λR, we have:
∫Tt#+ε{uFC,alt(t)C(Salt(t))+λS(t)S′alt(t)+λ′S(t)Salt(t)+λR(t)R′alt(t)+λ′R(t)Ralt(t)}dt−λS(T)Salt(T)+λS(t#+ε)Salt(t#+ε)−λR(T)Ralt(T)+λR(t#+ε)Ralt(t#+ε). |
Therefore, I2−I3 becomes
∫Tt#+ε{uFC(t)C(S(t))−uFC,alt(t)C(Salt(t))+λS(t)(S′(t)−S′alt(t))+λ′S(t)(S(t)−Salt(t))+λR(t)(R′(t)−R′alt(t))+λ′R(t)(R(t)−Ralt(t))}dt−λS(T)(S(T)−Salt(T))+λS(t#+ε)(S(t#+ε)−Salt(t#+ε))−λR(T)(R(T)−Ralt(T))+λR(t#+ε)(R(t#+ε)−Ralt(t#+ε)). | (B.3) |
Next, we will approximate the integrand of (B.3) with its Taylor series around the point (t,S(t),R(t),u(t)). We only keep track of the first-order terms of the expansion, because as ε→0, the integral already goes to zero. We have then that I2−I3 is approximately
∫Tt#+ε{(uFC(t)dCdS(S(t))+λS(t)dS′dS(t)+λR(t)dR′dS(t)+λ′S(t))(S(t)−Salt(t))∫Tt#+ε{+(λS(t)dS′dR(t)+λR(t)dR′dR(t)+λ′R(t))(R(t)−Ralt(t))∫Tt#+ε{+C(S(t))(uFC(t)−uFC,alt(t))+λS(t)dS′duSC(t)(uSC(t)−uSC,alt(t))∫Tt#+ε{+λR(t)dR′duRC(t)(uRC(t)−uRC,alt(t))+λS(t)dS′duSN(t)(uSN(t)−uSN,alt(t))∫Tt#+ε{+λR(t)dR′duRN(t)(uRN(t)−uRN,alt(t))}dt−λS(T)(S(T)−Salt(T))+λS(t#+ε)(S(t#+ε)−Salt(t#+ε))−λR(T)(R(T)−Ralt(T))+λR(t#+ε)(R(t#+ε)−Ralt(t#+ε)). | (B.4) |
Note that by the differential equations for λS and λR, and the fact that R(t#+ε)=Ralt(t#+ε), this simplifies to
∫Tt#+ε{C(S(t))(uFC(t)−uFC,alt(t))+λS(t)dS′duSC(t)(uSC(t)−uSC,alt(t))∫Tt#+ε{+λR(t)dR′duRC(t)(uRC(t)−uRC,alt(t))+λS(t)dS′duSN(t)(uSN(t)−uSN,alt(t))∫Tt#+ε{+λR(t)dR′duRN(t)(uRN(t)−uRN,alt(t))}dt+λS(t#+ε)(S(t#+ε)−Salt(t#+ε)). |
We can rewrite this in terms of the Hamiltonian:
I2−I3≈∫Tt#+ε{∂H∂uFC(t)(uFC(t)−uFC,alt(t))+∂H∂uSC(t)(uSC(t)−uSC,alt(t))∫Tt#+ε{+∂H∂uRC(t)(uRC(t)−uRC,alt(t))+∂H∂uSN(t)(uSN(t)−uSN,alt(t))∫Tt#+ε{+∂H∂uRN(t)(uRN(t)−uRN,alt(t))}dt+λS(t#+ε)(S(t#+ε)−Salt(t#+ε)). |
Now, because we know the structure of the optimal solution, namely that during the penultimate interval ∂H∂uFC(t)=∂H∂uSC(t),uSN=1,∂H∂uRC(t)=0=∂H∂uRN(t), and during the final interval uFC=1,∂H∂uFC(t)=C∗, we have that I2−I3 is approximately
∫t∗t#+ε{∂H∂uFC(uFC(t)+uSC(t)−uFC,alt(t)−uSC,alt(t))+∂H∂uSN(1−uSN,alt(t))}dt+∫Tt∗C∗(1−uFC,alt(t))dt+λS(t#+ε)(S(t#+ε)−Salt(t#+ε)). |
Assume that ε is chosen small enough that the optimal solution to the perturbed problem also begins in the penultimate interval at t=t#+ε. This gives us uFC,alt(t)+uSC,alt(t)=1=uSN,alt(t), and so
I2−I3≈C∗∫Tt∗uSC,alt(t)dt+λS(t#+ε)(S(t#+ε)−Salt(t#+ε)). |
Next, we use
∫t#+εt#F′(t)dt≈εF′(t#), |
which gives us
I1+I2−I3≈εF′(t#)+λS(t#+ε)(S(t#+ε)−Salt(t#+ε))+C∗∫Tt∗uSC,alt(t)dt. |
Using a linear approximation, and the fact that S(t#)=Salt(t#), we have that to O(ε)
λS(t#+ε)(S(t#+ε)−Salt(t#+ε))≈(λS(t#)+ελ′S(t#))(εS′(t#)−εS′alt(t#))≈ελS(t#)(S′(t#)−S′alt(t#)). |
Therefore, we have that
I1+I2−I3≈ε(F′(t#)+λS(t#)(S′(t#)−S′alt(t#))+C∗∫Tt∗uSC,alt(t)dt. | (B.5) |
The last step in obtaining the approximation in Theorem 1 is to show that the last term in (B.5) vanishes. To this end, we will show that t∗alt≤t∗, in which case the uSC,alt(t)=0 on [t∗,T]. To see why this is the case, first note that by the way we constructed the alternate growth path on [t#,t#+ε], we know that Salt(t#+ε)≥S(t#+ε). Because R(t)=R∗=Ralt(t) for all t in [t#,T], and Salt(t#+ε)>S(t#+ε), it must be the case that Salt(T)≥S(T), else there would exist a time t′ in (t#+ε,T) where S(t′)=Salt(t′) and R(t′)=Ralt(t′). In this case, the solutions to both the original problem and the perturbed problem would follow the same optimal trajectory for the remainder of the growing season. Now, because Salt(T)≥S(T) and by (3.11), we know that CS(t) is decreasing, we have that CS(Salt(T))≤CS(S(T)) which by (4.12) means that t∗alt≤t∗. Therefore, the last term in (B.5) vanishes and we are left with
I1+I2−I3≈ε(F′(t#)+λS(t#)(S′(t#)−S′alt(t#)). | (B.6) |
Now that we have obtained (B.6), it remains to show that I1+I2−I3>0, i.e., the optimal strategy given by the necessary conditions (shoot/fruit growth) is better than the strategy suggested by the marginal values (shoot-only growth). To this end, dropping the arguments after the first line,
I1+I2−I3≈ε(F′(t#)+λS(t#)(S′(t#)−S′alt(t#))=ε(uFCC+λS(S′−S′alt))=ε((1−uSC)C+λS(νSN∗G(uSCCνSN∗)−νSN∗G(CνSN∗))) |
Because here λSG′(uSCCνSN∗)=1, we have
I1+I2−I3≈ελS(CG′(uSCCνSN∗)−uSCCG′(uSCCνSN∗)+νSN∗G(uSCCνSN∗)−νSN∗G(CνSN∗))=ελSνSN∗(CνSN∗G′(uSCCνSN∗)−uSCCνSN∗G′(uSCCνSN∗)+G(uSCCνSN∗)−G(CνSN∗))=ελSνSN∗(CνSN∗G′(uSCCνSN∗)+G2(uSCCνSN∗)−G(CνSN∗)), |
where the last line follows from (3.20). At this point, showing that for a∈(0,1),x>0,
xG′(ax)+G2(ax)−G(x)>0 |
will be sufficient to show that I1+I2−I3>0. To that end, for a∈(0,1),x>0,
xG′(ax)+G2(ax)−G(x)=x(1+2ax)(1+ax+(ax)2)2+(ax)3(2+ax)(1+ax+(ax)2)2−x(1+x)1+x+x2=(a−1)2x3(1+2a+2ax+a2x)(1+ax+(ax)2)2(1+x+x2)>0. |
Therefore, we have the desired result, namely that
F(T)−Falt(T)≈ε(F′(t#)+λS(t#)(S′(t#)−S′alt(t#)))>0. | (4.22 revisited) |
In this appendix, we include the full proof of Theorem 2. Recall that t=ˆt denotes the transition point between the balanced growth phase and the penultimate interval, and z−SC and z+SC are the left and right limits of zSC at t=ˆt, respectively.
Theorem 2. zSC is continuous at t=ˆt.
Proof. First, note that by combining (5.26), (5.37), and (4.16), and writing ˆC=C(ˆt), we have that (x,y)=(z−SC,z+SC) must solve
C∗νSN∗G2(y)=(C∗−ˆC)(νSN∗G2(x)+ˆCG′(x)) | (5.50 revisited) |
C∗G′(y)=νSN∗G2(x)+ˆCG′(x). | (5.51 revisited) |
Because G2 and G′ are both one-to-one, we have that y is a function of x in (5.50) and (5.51). It will also be useful to equate the left-hand sides of (5.50) and (5.51) to obtain:
νSN∗G2(y)C∗−C=G′(y). | (C.1) |
We will break the argument up into three components.
First, we will show that for any given values of C∗,ˆC,νS, and N∗, equation (C.1) has only one positive real solution. To this end, setting k=νSN∗C∗−ˆC, equation (C.1) gives us
kG2(y)=G′(y)ky3(2+y)(1+y+y2)2=1+2y(1+y+y2)2ky3=1+2y2+y. |
Note that since C∗−ˆC≥0, we have that k>0. For any positive choice of k, it is a simple matter to see that the graphs of
x=ky3 and x=1+2y2+y |
have only one positive intersection point. Note that this tells us that G′(y)/G2(y) is invertible, and in particular, that there is only one admissible value of y (depending on k) that solves (C.1). Therefore, z+SC is completely determined by the values of C∗,ˆC,νS, and N∗.
Next, we note that x=y is a solution to equations (5.50) and (5.51). We will omit the details of this computation, as for both (5.50) and (5.51), it can be verified directly that x=y is a solution if and only if (C.1) holds, which we have already established.
Lastly, we will show that x=y is the only solution to (5.50) and (5.51). This, then, combined with the facts that (x,y)=(z−SC,z+SC) necessarily solves (5.50) and (5.51) and z+SC is completely determined by the values of C∗,ˆC,νS, and N∗ will establish that z−SC=z+SC, that is zSC is continuous at ˆt.
To show that x=y is the only solution to (5.50) and (5.51), we will use dydx for both these equations as well as the PCSU identity (3.21). We proceed with finding dydx for (5.50), and so as to avoid confusion, we will refer to y in this equation as yA. Differentiating (5.50) gives us
C∗νSN∗G′2(yA)dyAdx=(C∗−ˆC)(νSN∗G′2(x)+ˆCG″(x)) |
and so
dyAdx=(C∗−ˆC)(νSN∗G′2(x)+ˆCG″(x))C∗νSN∗G′2(yA)=−(C∗−ˆC)(ˆCG″(x)−νSN∗xG″(x))C∗νSyAG″(yA)=(C∗−ˆC)N∗G″(x)C∗yAG″(yA)(x−ˆCνSN∗). | (C.2) |
Now,
G″(x)=ddx1+2x(1+x+x2)2=−6x(1+x)(1+x+x2)3≤0 for x≥0 |
and so
(C∗−ˆC)N∗G″(x)C∗yAG″(yA)≥0 for x,yA≥0. |
In particular, this means that
sgn(dyAdx)=sgn(x−ˆCνSN∗). | (C.3) |
Next, we do the same thing for (5.51), here using yB to refer to y:
C∗G″(yB)dyBdx=νSN∗G′2(x)+ˆCG″(x) |
and so
dyBdx=νSN∗G′2(x)+ˆCG″(x)C∗G″(yB)=ˆCG″(x)−νSN∗xG″(x)C∗G″(yB)=−νSN∗G″(x)C∗G″(yB)(x−ˆCνSN∗). | (C.4) |
As before,
νSN∗G″(x)C∗G″(yB)≥0 for x,yB≥0, |
so, in particular, we have that
sgn(dyBdx)=−sgn(x−ˆCνSN∗) | (C.5) |
which means that
sgn(dyAdx)=−sgn(dyBdx). | (C.6) |
Additionally, note that both derivatives vanish at x=ˆCνSN∗. We will use these properties of dyAdx and dyBdx to show that x=y is the only admissible solution to (5.50) and (5.51).
Recall that (x,y)=(z−SC,z+SC) is necessarily a solution to (5.50) and (5.51), and we have established that the values of C∗,ˆC,νS, and N∗ determine a single value of y, and therefore z+SC, that solves both equations. Now, as z+SC=limt→ˆt+zSC, we have that
y=limt→ˆt+uSCCνSN∗≤ˆCνSN∗, |
which means that the y=x intersection point must be at or before x=ˆCνSN∗, the point where both dyAdx and dyBdx vanish. There are two possible cases to consider: either u+SC=1 or u+SC<1, where here u+SC=limt→ˆt+uSC(t).
1. Case 1: u+SC=1
In this case, the y=x intersection point occurs at (ˆCνSN∗,ˆCνSN∗). We have from (C.3) and (C.5) that x=ˆCνSN∗ is the unique global minimizer for yA and the unique global maximizer for yB, meaning that there are no additional intersection points.
2. Case 2: u+SC<1
In this case, the y=x intersection point occurs at (u+SCˆCνSN∗,u+SCˆCνSN∗), before x=ˆCνSN∗ where dyAdx=0=dyBdx. Now, we have from (C.3) and (C.5) that yA strictly decreases until x=ˆCνSN∗ and then strictly increases after that point, whereas yB strictly increases until x=ˆCνSN∗ and then strictly decreases after that point, implying the existence of exactly one more intersection point at some x=p>ˆCνSN∗. Now, as we established, there is only one solution for y, namely y=z+SC. Therefore, in this case, there is another intersection point (p,z+SC) where p>z+SC. Note that this is potentially admissible because we would have
p=z−SC=limt→ˆt−uSCCuSNνSN, |
which for admissible choices of the controls could be greater than u+SCˆCνSN∗. However, if z−SC=p, then as G′ is strictly decreasing, we have that G′(z−SC)<G′(z+SC), and so
λS(ˆt)G′(z−SC)<λS(ˆt)G′(z+SC)=1, |
which is impossible by (5.18). Therefore, this second intersection point is not admissible.
In both cases, we have shown that the only admissible intersection point for the curves defined by equations (5.50) and (5.51) is at the single intersection point where x=y, which because (x,y)=(z−SC,z+SC) is a solution means that zSC is continuous at ˆt.
This appendix includes a derivation, referenced in Section 6.5, for the differential equation for zSC and zRC during the balanced growth phase. In particular, we will derive equations (6.11) and (6.12), revisited below.
dzSCdt=NRνSνRG(zRC)G2(zSC)−CSG′(zSC)[νSG2(zSC)+zRCνRG′(zSC)]G″(zSC)[νSzSC−νRzRC] | (6.11 revisited) |
dzRCdt=NRνSG2(zSC)[G′(zRC)]2−CS[G′(zSC)]2G′(zRC)+G″(zSC)G′(zRC)dzSCdtG′(zSC)G″(zRC) | (6.12 revisited) |
In deriving these differential equations, we will use the notation zS=zSC and zR=zRC for simplicity. We begin by differentiating (5.18) and (5.19) with respect to t:
λ′SG′(zS)+λSG″(zS)dzSdt=λ′RG′(zR)+λRG″(zR)dzRdt | (D.1) |
λ′SνSG2(zS)+λSνSG′2(zS)dzSdt=λ′RνRG2(zR)+λRνRG′2(zR)dzRdt. | (D.2) |
Solving (D.1) for λRG″(zR)dzRdt gives us
λRG″(zR)dzRdt=λ′SG′(zS)+λSG″(zS)dzSdt−λ′RG′(zR). | (D.3) |
Likewise, making the substitution G′2(z)=−zG″(z) by means of (3.21), we can solve for zRνRλRG″(zR)dzRdt in (D.2) to get
zRνRλRG″(zR)dzRdt=λ′RνRG2(zR)−λ′SνSG2(zS)+zSλSνSG″(zS)dzSdt. | (D.4) |
Substituting (D.3) into (D.4) gives us
zRνR[λ′SG′(zS)+λSG″(zS)dzSdt−λ′RG′(zR)]=λ′RνRG2(zR)−λ′SνSG2(zS)+zSλSνSG″(zS)dzSdt. | (D.5) |
Solving for dzSdt, making use of (3.20), gives us
dzSdt=λ′S[νSG2(zS)+zRνRG′(zS)]−λ′RνRG(zR)λSG″(zS)[νSzS−νRzR]. | (D.6) |
Lastly, we use the differential equations for λS and λR during balanced growth, (5.22) and (5.23), to write both λ′S and λ′R in terms of λS. Canceling the resulting common factor of λS leaves us with
dzSdt=NRνSνRG(zR)G2(zS)−CSG′(zS)[νSG2(zS)+zRνRG′(zS)]G″(zS)[νSzS−νRzR] | (D.7) |
Similarly, we can use (5.22) and (5.23) along with (5.18) to simplify (D.3), which upon canceling the common factors of λS, gives us
dzRdt=NRνSG2(zS)[G′(zR)]2−CS[G′(zS)]2G′(zR)+G″(zS)G′(zR)dzSdtG′(zS)G″(zR). | (D.8) |
The authors declare that they have not used Generative-AI tools in the creation of this paper.
All authors declare no conflicts of interest in this paper.
This work was primarily conducted as part of David McMorris' PhD dissertation under the advisement of Glenn Ledder.
[1] |
Y. Iwasa, J. Roughgarden, Shoot/root balance of plants: Optimal growth of a system with many vegetative organs, Theor. Popul. Biol., 25 (1984), 78–105. https://doi.org/10.1016/0040-5809(84)90007-8 doi: 10.1016/0040-5809(84)90007-8
![]() |
[2] |
K. Velten, O. Richter, Optimal root/shoot-partitioning of carbohydrates in plants, Bull. Math. Biol., 57 (1995), 99–107. https://doi.org/10.1007/BF02458318 doi: 10.1007/BF02458318
![]() |
[3] |
B. J. Enquist, K. J. Niklas, Global allocation rules for patterns of biomass partitioning in seed plants, Science, 295 (2002), 1517–1520. https://doi.org/10.1126/science.1066360 doi: 10.1126/science.1066360
![]() |
[4] |
M. McCarthy, B. J. Enquist, Consistency between an allometric approach and optimal partitioning theory in global patterns of plant biomass allocation, Funct. Ecol., 21 (2007), 713–720. https://doi.org/10.1111/j.1365-2435.2007.01276.x doi: 10.1111/j.1365-2435.2007.01276.x
![]() |
[5] |
R. Dybzinski, C. Farrior, A. Wolf, P. B. Reich, S. W. Pacala, Evolutionarily stable strategy carbon allocation to foliage, wood, and fine roots in trees competing for light and nitrogen: An analytically tractable, individual-based model and quantitative comparisons to data, Am. Nat., 177 (2011), 153–166. https://doi.org/10.1086/657992 doi: 10.1086/657992
![]() |
[6] |
C. Feller, P. Favre, A. Janka, S. C. Zeeman, J.-P. Gabriel, D. Reinhardt, Mathematical modeling of the dynamics of shoot-root interactions and resource partitioning in plant growth, PLoS ONE, 10 (2015), e0127905. https://doi.org/10.1371/journal.pone.0127905 doi: 10.1371/journal.pone.0127905
![]() |
[7] |
G. Ledder, S. E. Russo, E. B. Muller, A. Peace, R. M. Nisbet, Local control of resource allocation is sufficient to model optimal dynamics in syntrophic systems, Theor. Ecol., 13 (2020), 481–501, https://doi.org/10.1007/s12080-020-00464-9 doi: 10.1007/s12080-020-00464-9
![]() |
[8] |
H. Poorter, K. J. Niklas, P. B. Reich, J. Oleksyn, P. Poot, L. Mommer, Biomass allocation to leaves, stems and roots: Meta-analyses of interspecific variation and environmental control, New Phytol., 193 (2012), 30–50. https://doi.org/10.1111/j.1469-8137.2011.03952.x doi: 10.1111/j.1469-8137.2011.03952.x
![]() |
[9] |
A. J. Bloom, F. S. Chapin III, H. A. Mooney, Resource limitation in plants-an economic analogy, Annu. Rev. Ecol. Evol. Syst., 16 (1985), 363–392. https://doi.org/10.1146/annurev.es.16.110185.002051 doi: 10.1146/annurev.es.16.110185.002051
![]() |
[10] |
H. Poorter, O. Nagel, The role of biomass allocation in the growth response of plants to different levels of light, co2, nutrients and water: A quantitative review, Aust. J. Plant Physiol., 27 (2000), 1191–1191. https://doi.org/10.1071/PP99173_CO doi: 10.1071/PP99173_CO
![]() |
[11] |
J. Weiner, Allocation, plasticity and allometry in plants, Pers. Plant Ecol., Evol. And Syst., 6 (2004), 207–215. https://doi.org/10.1078/1433-8319-00083 doi: 10.1078/1433-8319-00083
![]() |
[12] |
G. A. Fox, Annual plant life histories and the paradigm of resource allocation, Evolut. Ecol., 6 (1992), 482–499. https://doi.org/10.1007/BF02270693 doi: 10.1007/BF02270693
![]() |
[13] | S. Lenhart, J. T. Workman, Optimal control applied to biological models, Chapman and Hall/CRC, 2007. https://doi.org/10.1201/9781420011418 |
[14] | L. S. Pontryagin, E. F. Mishchenko, V. G. Boltyanskii, R. V. Gamkrelidze, The mathematical theory of optimal processes, Wiley, 1962. https://doi.org/10.1201/9780203749319 |
[15] | S. Kooijman, Dynamic energy budget theory for metabolic organisation, Cambridge university press, 2010. https://doi.org/10.1017/CBO9780511805400 |
[16] | M. I. Kamien, N. L. Schwartz, Dynamic optimization: the calculus of variations and optimal control in economics and management, North-Holland, 1991. |
[17] |
P. A. Hicks, Distribution of carbon/nitrogen ratio in the various organs of the wheat plant at different periods of its life history, New Phytol., 27 (1928), 108–116. https://doi.org/10.1111/j.1469-8137.1928.tb06735.x doi: 10.1111/j.1469-8137.1928.tb06735.x
![]() |
[18] |
V. Minden, M. Kleyer, Internal and external regulation of plant organ stoichiometry, Plant Biol., 16 (2014), 897–907. https://doi.org/10.1111/plb.12155 doi: 10.1111/plb.12155
![]() |
[19] |
S. Koshkin, Z. Zalles, M. F. Tobin, N. Toumbacaris, C. Spiess, Optimal allocation in annual plants with density-dependent fitness, Theory Biosci., 140 (2021), 177–196. https://doi.org/10.1007/s12064-021-00343-9 doi: 10.1007/s12064-021-00343-9
![]() |
[20] |
N. Chiariello, J. Roughgarden, Storage allocation in seasonal races of an annual plant: optimal versus actual allocation, Ecology, 65 (1984), 1290–1301. https://doi.org/10.2307/1938334 doi: 10.2307/1938334
![]() |
[21] | Y. Iwasa, Dynamic optimization of plant growth, Evol. Ecol. Res., 2, (2000), 437–455. |