Optimal seasonal schedule for the production of isoprene, a highly volatile biogenic VOC

The leaves of many trees emit volatile organic compounds (abbreviated as BVOCs), which protect them from various damages, such as herbivory, pathogens, and heat stress. For example, isoprene is highly volatile and is known to enhance the resistance to heat stress. In this study, we analyze the optimal seasonal schedule for producing isoprene in leaves to mitigate damage. We assume that photosynthetic rate, heat stress, and the stress-suppressing effect of isoprene may vary throughout the season. We seek the seasonal schedule of isoprene production that maximizes the total net photosynthesis using Pontryagin’s maximum principle. The isoprene production rate is determined by the changing balance between the cost and benefit of enhanced leaf protection over time. If heat stress peaks in midsummer, isoprene production can reach its highest levels during the summer. However, if a large portion of leaves is lost due to heat stress in a short period, the optimal schedule involves peaking isoprene production after the peak of heat stress. Both high photosynthetic rate and high isoprene volatility in midsummer make the peak of isoprene production in spring. These results can be clearly understood by distinguishing immediate impacts and the impacts of future expectations.

The amount of leaf area changes over time t following a differential equation.If the tree produces more BVOC on day t , the leaves are better protected from damages, resulting in improved photosynthesis from the surviving leaves.Whether this benefit exceeds the cost of producing the BVOC depends on various processes.Especially it depends on the day t .If t is close to the end of the growing season, the leaves do not have a high expectation of photosynthesis.and hence the value would be small.In contrast, in the early portion of the growing season, the same amount of leaf area would have a greater future prospect simply because the length of time from day t to the end of the season is longer.Hence, the isoprene production would be more profitable in the early portion of a season.
There are numerous processes influencing the benefit of isoprene production.To integrating them systematically, we employ a mathematical method developed in control theory, called Pontryagin's maximum principle 23,24 .This method has been adopted to elucidate the adaptive responses of organisms to environmental changes, explaining phenomena such as the timing of flowering within a year [25][26][27][28][29][30] , the balance of multiple organs of individual plants, such as shoot/root balance 31 , as well as in animal life history [32][33][34] .Additionally, it has been applied to growth of the dwarf males in barnacles 35 , production of shells for protecting mollusks from predation 36 , and the age-dependent concentration of defense chemicals, such as alkaloid, in biennial plants based on leaf age 37 .

Model
We consider a deciduous tree which produces a cohort of leaves synchronously on date t = 0 (leaf flushing).The leaves are gradually lost due to physical and biological damages.All remaining leaves are discarded at the end of the growing season ( t = T ).Let L(t) be the leaf area on day t .The dynamics of the leaf area are given as follows: for 0 < t < T .The initial value is L(0) = L 0 .The leaf area decreases due to random removal at rate u and the risk factor given by heat stress or herbivory h(t) .The intensity of the latter may vary over the season.The plant can reduce the risk by investing in BVOCs.We here consider a highly volatile compound, such as isoprene.Isoprene enhances the speed of removing harmful reactive oxygen species formed in leaves under the heat stress and helps the heat resistance of the leaves [5][6][7][8][9][10] .We consider the case in which the leaves do not maintain a large amount of isoprene for a prolonged period, assuming that all the isoprene produced on a day would be released within the same day and would not be carried over to the subsequent days, unlike terpenoids in Lamiacea or Pinacea species 38 .In Eq. (1), s(t) represents the amount of isoprene per leaf area produced on day t .In Eq. (1), the risk of leaf loss is reduced by factor 1/(1 + b(t)s(t)) , and the daily exponential rate of leaf area loss is u The objective function is to maximize the net photosynthesis of a cohort of leaves over the whole season, given as follows: Here p(t) represents the photosynthetic rate per unit leaf area, measured in terms of carbon gain.Equation ( 2) is the total amount of photosynthesis minus the investment in isoprene production.The amount of isoprene produced s(t) is measured in terms of the associated cost: the reduction of net CO 2 assimilation known to be caused by the higher allocation of leaf internal C sources (xylem sugars, starch, cytosolic sources) for isoprene biosynthesis 11 .This is a typical optimal control problem.The objective function is φ given in Eq. ( 2) and the control vari- able is s(t) for 0 < t < T , subject to the constraint 0 ≤ s(t) ≤ s max , where s max is the maximum rate of isoprene production per leaf area per day.Differential equation given by Eq. ( 1) is another constraint.We search for the optimal schedule of the isoprene production that maximizes φ.

Pontryagin's maximum principle
To apply Pontryagin's maximum principle to this optimal control problem, we define the Hamiltonian as follows: Here, (t) is a costate variable corresponding to the leaf area L(t) .Intuitively speaking, (t) implies the impact of a unit leaf area on day t in enhancing the net photosynthesis φ , as defined by Eq. (2) 24,31,39 .The costate variable follows the differential equation given by, We solve this equation in a backward manner, starting from the terminal condition specified as follows: The maximum principle dictates that the Hamiltonian must be maximized at each time t by selecting the control variable s(t) .Hamiltonian given in Eq. ( 3) is a function of s(t) with a negative second derivative.Accord- ing to the calculation in section A.1 of S1 Appendix A, we can derive that the optimal value maximizing the Hamiltonian given in Eq. ( 3) is equal to s(t) = s[ (t), t] , where s[ (t), t] is the quantity given as follows: (1) Equation ( 5) presents the optimal isoprene production rate, s(t) , as a function of (t).
To numerically obtain (t) , we can solve the differential equation given by Eq. ( 4) by integrating it from the end of the season, with the terminal condition (T) = 0 .At each time point t , we utilize s(t) = s[ (t), t] , as specified in Eq. ( 5).Refer to section A.2 of S1 Appendix A for a detailed derivation of the solution (t) including numerical procedure.
In the subsequent sections, we discuss how the optimal schedule of isoprene production depends on different parameters and rate functions included in the model.Understanding the results involves noting the dependence of isoprene production rate s(t) in the optimal schedule on various rate functions, as summarized in Fig. 1.
First, Eq. ( 5) indicates that the isoprene production rate s(t) is determined by strength of heat stress h(t) , the effectiveness of isoprene in suppressing stress b(t) , and the marginal value of leaf area (t) on the same day t .This implies that the changes in these quantities affect the isoprene production rate immediately, which may be called "immediate impacts." In Fig. 1, this is indicated by three red arrows from h(t) , b(t) , and (t) to s(t) , respectively.
Second, (t) is obtained from integrating differential equation given by Eq. ( 4) from t to the end of the season T .Hence, (t) depends on the photosynthetic rate, heat stress, effectiveness of BOVC after t (i.e., p(t′) , h(t′) and b(t′) for t < t′ < T ).For example, an enhancement of photosynthestic rate p(t′) on day t′ affects the value of leaf area (t) on an earlier date t (i.e.,t < t′ ).This indicates that the value of leaf area (t) depends on the expected value of photosynthetic rate in the future.We may call this as "the impact of future expectations." In Fig. 1, this is indicated by a blue arrow from p(t) to (t).
Photosynthetic rate p(t) does not affect s(t) directly, but it affects s(t) indirectly via modifying marginal value (t) .In Fig. 1, this is shown by a blue arrow from p(t) to (t) together with a red arrow from (t) to s(t) .Both h(t) and b(t) exert direct and indirect effects (i.e., immediate impacts and the impacts of future expectations) on the optimal isoprene production s(t).
In section A.3 of S1 Appendix A we derive the differential equation Eq. ( 4) only from the definition of (t) as the value of unit leaf area on day t .We also justify the argument that the optimal value of isoprene production rate is given by Eq. ( 5) without resorting the official mathematical formalism of the Hamiltonian.

Optimal seasonal schedule
Based on the method explained above, we generated the optimal seasonal schedule.

Final phase without investing isoprene production
As a season progresses, the marginal value of leaf area (t) decreases following Eq. ( 5 and reaches zero at t = T .Toward the end of the growing season, there is a period during which no isoprene is produced.This is because protecting leaves from damage is not profitable when the number of remaining days is small.Let t s be the date at Figure 1.Scheme of how the optimal isoprene production is affected by time-dependent rates.The optimal isoprene production rate s(t) is determined by heat stress h(t) , the effectiveness of BVOC b(t) , and the marginal value of leaf area (t) on the same day.These immediate impacts are indicated by red arrows.Additionally, s(t) is influenced by future values of p(t), h(t) , and b(t) , which modify (t) on day t .These impacts of future expectations are indicated by blue arrows.In the text, we refer to these influences as direct and indirect effects, respectively.
Vol:.( 1234567890) www.nature.com/scientificreports/which investment to isoprene production ceases.We integrate Eq. ( 4) with s(t) replaced by s[ (t), t] in a backward manner from the final date T with (T) = 0 .From (t s )b(t s )h(t s ) = 1 , we have the following condition for t s .
When photosynthesis rate is slow, the isoprene is not effective in suppressing stress, or heat stress is small, the isoprene production is not profitable to the tree.Equation ( 6) is not satisfied, because the left-hand side of Eq. ( 6) is larger than the right-hand side for all t s .When all the rates are independent of time, Eq. ( 6) specifies the onset of the final phase of no isoprene production as follows: The derivation is explained in S1 Appendix A. Equation (7) indicates that the isoprene production lasts longer if the photosynthesis is fast (large p ), heat stress is strong (large h ), and isoprene is effective (large b).

When all the rates are independent of time t
First, we consider simple situations in which all rates (photosynthesis p(t) , harm h(t) , and effectiveness of isoprene b(t) ) are independent of time, t .Figure 2 illustrates a typical case in which there are three phases.In the first phase, just after leaf flush, isoprene production occurs at the fastest rate.This period is followed by the second phase, during which isoprene is produced at an intermediate rate (between 0 and s max ).After t s , there is a final period in which no isoprene is produced.The switching time t s is given by Eq. ( 7).
The level of isoprene production is also influenced by parameters within the model.For example, while the seasonal pattern of the isoprene production s(t) remains consistent, an increase in the photosynthetic rate results in a higher s(t) value and maintains positivity for an extended duration within the season (refer to Fig. S1 in SI Appendix B for detail).In SI Appendix B, we explain the results of parameter dependence in detail.When s max is sufficiently large, the optimal schedule does not include Phase I, indicating that isoprene production is not constrained by the maximum daily production rate.Then, the maximum rate of isoprene production is s(0) .The elasticities of s(0) with respect to p , h , b , and T were positive and of a similar order of magnitude, although their precise values vary depending on the choice of the standard parameter sets.In contrast, the elasticity with respect to u was negative and significantly smaller in magnitude than the others.From these results, we concluded that the optimal isoprene production tends to increase with a higher photosynthesis rate, increased heat stress, improved efficiency of the chemical, and a longer growing season.

When harm h(t) peaks in the summer
Heat stress is expected to be more pronounced in the summer than in spring or fall.Figure 3 illustrates the cases where h(t) reaches its peak in the middle of the growing season.Isoprene production is low or small in the spring or early summer, increases to a peak level in the mid-summer, and decreases in the late summer.The heat stress is an exponentially modulated cosine function, h(t) = exp[a − bcos{2πt/T}] , which attains the maximum exp[a + b] at t = T/2 and the minimum exp[a − b] at t = 0 and t = T .A larger b indicates a sharper peak of heat stress.We may call b as "shape factor." (i) High isoprene production during the period of high heat stress Optimal schedule when all parameters are constant.The BVOC is produced at the maximum rate s(t) = s max in the first period, an intermediate rate 0 < s(t) < s max is the second period, and no BVOC is produced in the final period s(t) = 0 .Parameters are: Figure 4 shows a numerical example of the optimal schedule of isoprene production, exhibiting three different values of strength of heat stress.The optimal isoprene production had a peak in the mid-summer.This is plausible because isoprene functions to cope with heat stress in the mid-summer.The three parts are for cases with the same shape constant b but different a .As a increases, the total heat stress t H total = T 0 h(t)dt increases, and the peak date of isoprene production remained near the peak date of heat stress (midsummer).(ii) Post-risk enhancement However, if s(t) curve is closely examined, the date of the highest rate of isoprene production occurred a little later than the peak of the heat stress.The difference between the two peaks became more pronounced for a larger total heat stress (Fig. 4).This suggests that the isoprene is produced at the maximum rate in late summer or even in early autumn.The reason is explained by the non-monotonic shape of the marginal value of leaf area (t) .When the heat stress was very strong, (t) became the lowest before the mid-summer of heat stress and increased rapidly in the high-risk period, resulting in a rather high value of leaf area in the late summer or autumn, as shown in Fig. 4. The marginal value of leaf area exhibited this behavior because leaf area sampled before the high-risk period was more likely to be removed during that time, rendering the leaf area sampled after the risky period more value than that sampled before it.We call this increase of the value as "post-risk enhancement."The magnitude of this effect was related to the amount of leaf areas loss in this peak stress, observed in the decrease in L(t).In SI Appendix C, we explain the results for different constant b for the shape of h(t) , with the total heat risk H total = T 0 h(t)dt controlled, which suggest that the magnitude of the peak shift was more strongly controlled by total heat risk H total than by shape constant b.

Effects of seasonal changes in photosynthesis and volatility
We also examined the effect of seasonal change in photosynthetic rate of leaves p(t) and the effectiveness of isoprene in suppressing heat stress b(t).
Figure 5(a) illustrates the results when p(t) has a peak in midsummer, with other rates constant.The opti- mal isoprene production reaches its peak in spring or in early summer and decreases sharply in midsummer and autumn.This result can be understood from the pattern of marginal value of leaf area (t) , which declines  .All other parameters are constant.We fixed a 1 = 0 .Three cases differ in b 1 , which changes H total ( b 1 = 1 , 2 , and 4 , respectively in top, middle, and bottom parts).As a increases, the total stress increases ( H total = 1.27 , 2.28 , and 11.3 , respectively).Isoprene production peaks slightly later than the peak heat stress, with the magnitude of their difference increasing with the total stress.Other parameters are: p = 1 , b = 2 , u = 0.0001 , T = 1.rapidly in summer.As shown in Fig. 1, photosynthetic rate p(t) affects the isoprene production s(t) indirectly via modifying the marginal value (t) .A high peak p(t) in midsummer leads to a high peak of (t) in spring and the decrease of (t) and s(t) in summer, because p(t) has only the impact of future expectations on s(t) (see Fig. 1).
In midsummer, the volatility of the chemicals may increase due to high temperatures, leading to a drop of b(t) in midsummer.Figure 5(b) illustrates the effect of this seasonality.The optimal isoprene production is high in spring or in early summer, but it becomes zero in midsummer and in autumn.b(t) has both direct and indirect effects on s(t) (see Fig. 1).The drop of b(t) in midsummer would directly reduce s(t) in summer (immediate impacts).Furthermore, the indirect effect of b(t) drop causes a rapid decline in (t) (impacts of future expecta- tions).The combination of these effects may stop isoprene production in midsummer and autumn, leading to the peak of s(t) in spring.

Discussion
Many trees produce and emit biogenic volatile organic compounds (BVOCs), which protect their leaves from various damages, including herbivory, pathogens, and heat stress 1 .The amount of isoprene production by trees is very large 2 , and BVOCs may modify the local and global climates 3,4,40 .Studies of the molecular and physiological mechanisms have also been advanced 21,[41][42][43][44][45][46] .Understanding the seasonal pattern of gene expression would be helpful in revealing the mechanism controlling genes, considering the respective roles of their products.
In this study, to understand the adaptation of different seasonal patterns of isoprene production, we analyzed the optimal seasonal schedule for producing highly volatile BVOCs, such as isoprene, in leaves that mitigate damage caused by heat stress [5][6][7][8][9]13 . We ormalized a dynamic optimization problem for the seasonal schedule of Isoprene production.Using Pontryagin's maximum principle, we searched for the optimal seasonal schedule that maximizes the total net photosynthesis and examined how it varies with seasonal patterns of photosynthetic rate, heat stress, and the stress-suppressing effect of isoprene.

Isoprene production level depends on future expectations
Producing isoprene provides benefits to the tree through increased resistance to heat stress and elevated defense levels against herbivores/pathogens. Since the benefit arises from the preservation of leaves, it depends on the length of the period for them to perform photosynthesis.Toward the end of the growing season, the tree makes no BVOC because the benefit is smaller than the cost.If the time until the end of the season is sufficiently long, producing some isoprene becomes beneficial.To represent how the value of preservation of leaves changes with time in the growing season, marginal value of leaf area (t) is defined.In section A.4 of SI Appendix A, we derive that (t) is equal to the total net production to be made until the end of the season per leaf area.
The second aspect important in evaluating the benefit of additional isoprene production is that the effect of stress suppression by additional amount diminishes as stress is reduced by the chemical already produced.The benefit of isoprene production of an additional unit amount is the effect of suppressing leaf area loss multiplied by the marginal value of leaf area.Maintaining a production rate that strikes the balance between the benefits and costs of isoprene production provides the optimal solution.
We may ask the amount and seasonal pattern of isoprene production that is the most profitable to the tree.In simple cases where the photosynthetic rate, stress intensity, and the stress-mitigating effect of isoprene are all constant, isoprene production should be high in the early growing season, decreasing gradually and reaching zero toward the end of the season, as shown in Fig. 3.According to elasticity analysis, trees tend to produce more isoprene in highly productive environments, under increased heat stress, with more effective inhibition by isoprene, during a longer growing season, and with a smaller leaf loss rate.

Seasonal changes of different processes shape the isoprene seasonal pattern
Next, we investigated the optimal isoprene production schedule in under seasonal changes.When heat stress h(t) peaks in the mid-summer with others constant, optimal isoprene production has a peak in midsummer, as shown in Figs. 4, S1.This result is plausible.
The bottom parts of Fig. 5 illustrated the results when the total amount of stress T 0 h(t)dt was large.The maximum rate of isoprene production per unit leaf area occurs later than the peak of heat stress h(t) , namely, late summer or autumn, rather than midsummer.This counter-intuitive result is explained by the nonmonotonic seasonal change in the marginal value of leaf area: (t) is small just before the peak heat stress.It rises sharply during the mid-summer and has the maximum after the peak of h(t) .The direct effect of seasonal change of h(t) makes s(t) high in mid-summer, while the indirect effect of h(t) causes the delay in the peak s(t) .The seasonal pattern of s(t) is the combined result of these two effects (see Fig. 1).We called the increase in (t) during the time of a strong heat stress "post-risk enhancement" of the marginal value of leaf area.This effect is pronounced only when a substantial fraction of leaves is lost during the period.While having peak activity in late summer seems counter intuitive, but it is interesting to note that Mayrhofer et al. 41 reported in grey popular leaves that the protein concentration and enzymatic activity of Isoprene synthase (PcISPS) peaked in late summer, while transcriptional levels of PcISPS gene were the highest in early summer.
When photosynthetic rate p(t) peaks in midsummer with other rates constant, the optimal isoprene produc- tion has a peak in spring.The summer peak of p(t) makes the marginal value of leaf area (t) decline quickly in midsummer and becomes small in autumn.Since p(t) has no immediate impacts on s(t) (see Fig. 1), s(t) has a peak in spring.In a similar manner, when the effectiveness of the isoprene in mitigating heat stress b(t) drops in midsummer due to the high volatility, the optimal isoprene production has a peak in spring.These are not the pattern observed for genes involved in the pathway for production of isoprene 21 .Hence, we may conclude that the simplest way to explain the peak isoprene production in midsummer is the high heat stress.The seasonal changes in photosynthetic rate or volatility play a less significant role.To comprehend how the seasonality of different processes shapes the optimal isoprene production pattern, we need to distinguish between the immediate impacts and those stemming from future expectations (refer to Fig. 1).

Future problems
There are several different ways of future theoretical studies extending the model examined in this paper.
(i) mildly volatile BVOC In a study of genome-wide gene expression data, focusing on Quercus glauca and Lithocarpus edulis in warm temperate regions of Japanese islands, genes related to sesquiterpene production were highly expressed in spring, while genes related to isoprene production exhibited high expression in the middle of the summer 21,22 .Sesquiterpene is much less volatile than isoprene.The observed difference in the seasonal pattern of gene expression could be attributed to the difference in the turnover rate of the respective chemicals.To address this hypothesis, we may consider the optimal seasonal schedule of producing BVOC with mild (or low) volatility.We conjecture that the most efficient seasonal schedule could involve concentrating the production of most chemicals at the beginning of the growing season, without additional production in the rest of the year.To be specific, we may assume that the risk is reduced by the presence of BVOC per leaf area V (t) , as h(t)/(1 + aV (t)) .The amount of BVOC on day t is not proportional to the production on the same day but it is determined by the amount produced on or before the day.For example, V (t) might follow a differential equation, dV /dt = s(t) − k(t)V (t) , where s(t) is the production rate of BVOC.The model is more challeng- ing to analyze than the one discussed in the present study due to the presence of two differential equations in the system.It is a problem worthy of being addressed in the future studies.(ii) BVOC emission as communication between and within individuals Plants under stress might benefit from signaling to other individuals.For instance, when one tree is attacked by herbivores, it may release monoterpenes and sesquiterpenes and attract predatory insects, which help in eliminating the herbivorous insects [14][15][16][17][18] .Additionally, when a tree is under attack, it may emit defence-related BVOCs, makes the recipient trees start producing the BVOCs.This could induce a form of defense in nearby trees, known as "induced systemic resistance" 47 .Furthermore, when a branch is attacked by herbivores and emits BVOCs, other branches of the same individual initiate the production of the BVOCs, thereby enhancing the overall the defense level.The diverse aspects of communication, among different species, among individuals of the same species, and within a single individual, will be an important subject of study in theoretical biology.(iii) Alternative method of calculating optimal control In the present study, Pontryagin's maximum principle was adopted for analysis.Dynamic programming is an alternative method of solving the dynamic optimization problem 48,49 , which has been used in the evolutionary ecology to analyze the optimal life history schedule of plants and animals 39,[50][51][52] .Since dynamic programming easily extends to situations including discrete-time dynamics and stochastic environments, developing corresponding analyses based on the dynamic programming method may provide us with a different perspective on the problem.(iv) Mild heat stress reducing the photosynthetic rate in a reversible manner A severe heat stress is known to damage the photosynthetic machinery in an irreversible manner 12 , which can be effectively modeled as the loss of leaf area as in this paper.However, mild heat stress is known to reduce the photosynthetic rate in a reversible manner 12 .To model such situations, we need to consider the possibility that the photosynthetic rate can be reduced by the heat, which may also be mitigated by isoprene production.For example, p(t) in Eq. ( 2) is replaced by p(t)exp[−γ h(t)/(1 + b(t)s(t))] .In this altered situation, the optimal isoprene production s(t) may depend on p(t) .This is an important theme for future studies.The model studied in this paper was chosen to be the simplest.We may extend it in several different directions, incorporating various complexities, such as unpredictable environmental fluctuations like irregular heatwaves, pest outbreaks, or weather patterns, and interactions with a broader ecosystem.Additionally, we will consider the resilience and adaptability of the BVOC production strategy and the potential feedback loops between BVOC emissions and climate change.This includes examining both the local cooling effects (due to aerosol formation) and the global warming potential (due to methane and ozone formation).

Figure 4 .
Figure 4. Optimal schedule of isoprene production when heat stress has a peak in the middle of the season (i.e.summer).We adopted an exponentially modulated cosine function, h(t) = exp a 1 − b 1 cos 2πt T