Towards a unified theory of plant photosynthesis and hydraulics

The global carbon and water cycles are governed by the coupling of CO2 and water vapour exchanges through the leaves of terrestrial plants, controlled by plant adaptations to balance carbon gains and hydraulic risks. We introduce a trait-based optimality theory that unifies the treatment of stomatal responses and biochemical acclimation of plants to environments changing on multiple timescales. Tested with experimental data from 18 species, our model successfully predicts the simultaneous decline in carbon assimilation rate, stomatal conductance and photosynthetic capacity during progressive soil drought. It also correctly predicts the dependencies of gas exchange on atmospheric vapour pressure deficit, temperature and CO2. Model predictions are also consistent with widely observed empirical patterns, such as the distribution of hydraulic strategies. Our unified theory opens new avenues for reliably modelling the interactive effects of drying soil and rising atmospheric CO2 on global photosynthesis and transpiration. Using trait-based optimality theory that unifies stomatal responses and acclimation of plants to changing environments, this study builds a model of the coupling of CO2 and water vapour exchanges through the leaves. This successfully predicts the simultaneous decline in carbon assimilation, stomatal conductance and photosynthetic capacity during progressive droughts.

The global carbon and water cycles are governed by the coupling of CO 2 and water vapour exchanges through the leaves of terrestrial plants, controlled by plant adaptations to balance carbon gains and hydraulic risks. We introduce a trait-based optimality theory that unifies the treatment of stomatal responses and biochemical acclimation of plants to environments changing on multiple timescales. Tested with experimental data from 18 species, our model successfully predicts the simultaneous decline in carbon assimilation rate, stomatal conductance and photosynthetic capacity during progressive soil drought. It also correctly predicts the dependencies of gas exchange on atmospheric vapour pressure deficit, temperature and CO 2 . Model predictions are also consistent with widely observed empirical patterns, such as the distribution of hydraulic strategies. Our unified theory opens new avenues for reliably modelling the interactive effects of drying soil and rising atmospheric CO 2 on global photosynthesis and transpiration.
The fundamental dilemma of plants following the C3 photosynthetic pathway is that when stomata, that is, the tiny 'valves' on the surface of leaves, are opened to take in carbon dioxide (CO 2 ) for carbon assimilation, water is lost through them via transpiration 1 . The plant's transpiration stream is maintained by negative water potentials (suction pressures) in roots, transport vessels and leaves. Withstanding negative water potentials requires adapted stem, leaf and root tissues or energy-intensive repair efforts, and extreme water potentials in the xylem can lead to hydraulic failure [2][3][4] . The risks of hydraulic failure increase when water availability declines across the plants' rooting zone or when vapour pressure deficit increases at their leaf surfaces. Plants can avoid hydraulic failure by closing their stomatal openings in response to dry soil and atmospheric conditions. However, closing the stomata also leads to a decline in carbon assimilation, creating a tight coupling between carbon uptake and water loss. At the ecosystem level, this coupling of the carbon and water cycles governs the rates of gross primary production (GPP) and evapotranspiration in response to water stress. On one hand, rising atmospheric CO 2 and increased precipitation are enhancing water use efficiency 5,6 , potentially increasing tree growth rates. On the other hand, rising atmospheric vapour pressure deficits are leading to decreases in stomatal conductance 7 , and rising frequency and intensity of droughts are leading to increased mortality rates 8 . It has been argued that a persistent increase in tree mortality rates, together with a saturating increase in growth rates, is negatively affecting the carbon sink of tropical forests 9 . Accurate predictions of carbon and water fluxes under water stress thus require vegetation models that explicitly account for plant hydraulic processes 10 Article https://doi.org/10.1038/s41477-022-01244-5 to resolve the limiting effect of atmospheric water demand and soil moisture stress on plant photosynthesis 11 .
A plant's hydraulic machinery places key constraints on how much water it can transpire and, consequently, on its stomatal conductance. Considerable effort has gone into the development of stomatal control models with an explicit treatment of plant hydraulics (see reviews 12,13 ). Hydraulically explicit stomatal models have shown success in simulating short-term stomatal responses to drying soil and air on sub-daily and daily timescales [14][15][16][17] and are now being implemented in Earth System Models [18][19][20][21] . However, we still lack understanding of how plant physiology acclimates to the development of soil-moisture drought on daily to weekly timescales and how such longer-term acclimation in turn affects stomatal sensitivity to short-term water stress. Such understanding is especially crucial for predicting stomatal and biochemical responses to novel environments and for explaining widely observed patterns related to plant hydraulic strategies (Box 1) in a parsimonious way.
The classic stomatal optimization model 22 states that plants adjust their stomatal conductance to maximize total carbon assimilation for a fixed amount of water loss, by assuming a constant unit cost for transpired water. This model implies that plants can save water for future use. However, recent stomatal models recognize that plants competitively consume available water 23 . Therefore, an alternative approach conceives the costs of transpiration as arising from the risks of hydraulic failure and the structural and energetic expenditures for withstanding high suction pressures. Thus, many extensions of this classic model explicitly represent plant hydraulics and the associated costs [23][24][25] . These models require an a priori specification of photosynthetic capacity, which then becomes an additional parameter to be fitted to enable accurate predictions of assimilation rates. By contrast, the least-costs optimization framework 26 includes the costs of maintaining carboxylation capacity, reflecting a trade-off between investing in photosynthetic and hydraulic capacities 27 . Building upon this approach, a recent model predicts acclimated carboxylation capacity 28 using the photosynthetic-coordination hypothesis 28,29 . It also explicitly optimizes electron-transport capacity (albeit using a separate optimization criterion) 28 , and has been successful in predicting CO 2 assimilation rates and leaf-internal CO 2 concentrations across climatic gradients. However, this model requires an empirical factor to account for the effects of soil moisture [30][31][32] .
Here we develop a unified first-principles theory combining the photosynthetic-coordination hypothesis with the principles of plant hydraulics within a single optimality framework. Our framework simultaneously predicts the stomatal responses and biochemical acclimation of leaves to environments changing on multiple timescales. We test the resulting model predictions with published data obtained from soil drought experiments conducted with 18 plant species spanning diverse plant functional types. We show that, with just three hydraulic traits and two parameters, our model correctly predicts key observations related to plant photosynthetic responses and hydraulic strategies, as described in Box 1.

Model summary
We now list the principles and hypotheses underlying our model in general terms, followed by a summary of the optimality framework, plant traits used in our model, the interpretation of model parameters and our strategy for testing the model with experimental data. A detailed model description is presented in Methods, and a full derivation of the model is presented in Supplementary Information section 1.

Model principles and hypotheses
Our model is based on three principles and hypotheses as follows.
(1) Water-balance principle. Any plant must maintain a continuous stream of water across its entire hydraulic pathway (through roots, stems and leaves) to ensure that the atmospheric demand

Box 1
Widely observed empirical patterns in leaf photosynthetic responses and plant hydraulic strategies as targets for model-based predictions Stomatal and biochemical responses to soil and atmospheric drought 1. As soil moisture decreases or vapour pressure deficit (D) increases, the first response of leaves is to reduce their stomatal opening to alleviate water stress. Since carbon uptake and water loss occur through stomata, photosynthesis and transpiration both decline with stomatal closure and thus, with decreasing soil moisture 30 . 2. As assimilation declines, maintaining photosynthetic capacity becomes increasingly unprofitable. Therefore, in the short term, leaf photosynthetic capacity also declines with decreasing soil moisture 44,45 . Yet, in the long term, plants acclimate by shedding their leaves, which reduces transpiration demand and allows assimilation to recover. Thus, in the long term, high photosynthetic capacity can be maintained at the expense of a reduced foliage surface area 58 . 3. As D increases, the leaf-internal-to-external CO 2 ratio (χ) declines. Various functional forms have been used to describe this decline, these forms having been mostly derived from limited empirical data 77 . A widely used relationship, predicted by simple stomatal optimization models 26,78 , is χ = ξ/(ξ + √D), where D is expressed as a fraction of atmospheric pressure and ξ is a constant. This implies that logit(χ) varies linearly with log(D) with a slope of −0.5, a value often targeted by modellers 23,28 . However, a recent study 46 analysing data from hundreds of species along aridity gradients has reported slope values of −0.76 ± 0.15, with remarkable consistency across species. The extent and consistency of these observations suggest that this value ought to be taken as a new target for model-based predictions.
Hydraulic strategies and trait-adaptations 4. As soil moisture decreases or D increases, xylem water potentials become increasingly negative. Extremely negative water potentials create embolisms in the xylem, which have been linked to increased risks of plant mortality ('hydraulic failure'). To avoid these risks, plants close their stomata before the onset of substantial xylem embolism 4,37,49,50,79 67,80 . At one end are isohydric (drought-avoiding) species that maintain a constant leaf water potential by closing the stomata as soil water potential decreases, at the cost of reduced carbon assimilation. At the other end are extreme anisohydric (drought-tolerating) species that keep their stomata open even in the face of decreasing soil water potential to maintain high CO 2 uptake, at the risk of hydraulic failure. In between are isohydrodynamic species that maintain a relatively constant soil-to-leaf water-potential difference.
Article https://doi.org/10.1038/s41477-022-01244-5 for transpiration is met by the supply of water from the soil 33 . If supply does not equal demand, xylem may embolize or leaves and roots may get damaged, causing catastrophic failure of the hydraulic system. Demand through transpiration depends on the stomatal conductance g s and the atmospheric vapour pressure deficit D, whereas supply depends on the soil-to-leaf water-potential difference Δψ and the hydraulic properties of the transpiration pathway. Therefore, this principle predicts g s as a function of Δψ and is widely used in stomatal models that explicitly represent water transport. We use the term 'principle' rather than 'hypothesis' for this assumption to indicate its rooting in basic physical laws. (2) Photosynthetic-coordination hypothesis. Photosynthetic carbon assimilation is limited by a plant's capacity for carboxylation V cmax and by light availability I abs , which, together with the electron-transport capacity J max , determine the rates of biochemical and photochemical reactions governing CO 2 fixation 34 . In general, the rate of photosynthesis is the minimum of the carboxylation-limited rate A c and the light-limited rate A j . The light-limited rate is further modulated by J max . Since the carboxylation and electron-transport capacities are costly to maintain, they are hypothesized to acclimate to typical daytime conditions on a weekly timescale, such that the two photosynthetic rates are coordinated, that is, A c = A j 29,35 . (3) Profit-maximization hypothesis. We posit that, on a weekly timescale (medium-term responses), plants simultaneously optimize their photosynthetic capacity and stomatal conductance to maximize net assimilation (profit, F), after accounting for the costs of maintaining photosynthetic capacity and the hydraulic pathway, including the risks of hydraulic failure. On a daily timescale (short-term responses), the acclimated photosynthetic capacities are fixed, and plants can optimize only their stomatal conductance. The parameters scaling the photosynthetic and hydraulic costs, α and γ, respectively, are the only two latent (that is, not directly observable) parameters in our model and are henceforth called 'unit costs'.

Hydraulic pathway and hydraulic traits
Water from the soil first enters the roots, where it flows through the root cortex, the endodermis and the stele. It then flows through the xylem in roots, the stem and leaf veins. After exiting the xylem, it flows through the bundle sheath and spongy mesophyll cells in the leaf, until it evaporates from the stomatal cell walls and diffuses out to the ambient air 36 (Fig. 1a). Under moderate water stress, which is the focus of the present work, the outside-xylem segments of the pathway may experience reversible losses in conductivity. For example, the soil-root interface may lose conductivity as roots shrink and disconnect from the soil, while leaves and roots may lose conductivity due to reductions in aquaporin activity or cell membrane permeability 37,38 . Extremely high suction pressures may lead to irreversible loss of xylem conductivity due to cavitation, although some species can reverse it by refilling xylem conduits or growing new xylem 39 . The flow of water through the plant can be described as depending on three effective hydraulic traits, which characterize the combined effect of the individual segments of the pathway: (1) the maximum plant conductance per unit leaf area, that is, the maximum leaf-specific whole-plant conductance K p , (2) the water potential ψ 50 that causes a 50% loss of whole-plant conductance, and (3) a shape parameter b that determines the sensitivity of conductance loss to water potential during progressive drought. There is increasing evidence that roots are the most resistive parts of the hydraulic pathway 40 and leaves  Fig. 1 | Schematic representation of our model, underlying first principles and notation. a, Water-transport pathway. Purple labels indicate the three hydraulic traits that determine the conductance to water flow of each of the three segments of the water-transport pathway. Water potentials are shown at various points along the pathway: ψ s in soil, ψ r in roots at the beginning of the xylem segment, ψ x at the end of the xylem segment and ψ l in leaves near the stomata. The soil-to-leaf water potential difference Δψ = ψ s − ψ l thus comprises the successive pressure drops along the three segments, that is, Δψ r = ψ s − ψ r along the radial outside-xylem segment within the roots, Δψ x = ψ r − ψ x along the xylem and Δψ l = ψ x − ψ l along the outside-xylem segment within the leaves. b, Model-calibration pathway. The model takes as inputs three effective whole-plant hydraulic traits (K p , ψ 50 and b) together with two cost parameters (the unit costs of photosynthetic and hydraulic capacities, α and γ, respectively). It predicts as outputs the optimal values (denoted by asterisks) of stomatal conductance g * s , assimilation rate A * , transpiration E * , acclimated photosynthetic capacities V * cmax and J * max , soil-to-leaf water-potential difference Δψ * and leaf internal-to-external CO 2 ratio χ * . Each variable is first calculated as a function of Δψ and χ, as shown by the four light-green arrows, from which the optimal combination (Δψ * , χ * ) is then calculated by maximizing profit F according to equation (1). Blue arrows and boxes indicate the process through which the best-fit traits and unit costs for each species are calculated by minimizing the model error. Orange labels indicate the three principles and hypotheses underlying the model, displayed next to the processes they affect.
Article https://doi.org/10.1038/s41477-022-01244-5 are the most vulnerable 37,41 . This is primarily due to the properties of the outside-xylem segments in roots and leaves, which thus form the hydraulic bottleneck of a plant. Further, there are lags in the recovery of root and stem conductivity upon rehydration, causing hysteresis in the response of the conductivity of these tissues to water potential 42 . While an explicit treatment of such hysteresis may be important for predicting root recovery on short timescales and xylem recovery from extreme drought, we focus on one-way drying in this study, reflecting data limitations and maintaining simplicity and analytical tractability.

Medium-term responses
To predict the acclimation of photosynthetic capacity on a weekly timescale, we assume that plants independently control their weekly-average stomatal conductance g s and their electron-transport capacity J max to maximize their net profit F, as defined below. After expressing all quantities in Fig. 1b in dependence on g s and J max , or equivalently, in a mathematically more convenient form, in terms of the leaf internal-to-external CO 2 ratio χ and the soil-leaf water-potential difference Δψ (Methods and Supplementary Information section 1.3), F can be written as where A is the assimilation rate calculated by combining the standard biochemical model of photosynthesis 34 with the photosynthetic-coordination hypothesis (equation 6 in Methods). We find the optimal solution (χ * , Δψ * ) semi-analytically by first calculating the derivatives of F with respect to χ and Δψ analytically (Supplementary equation 16) and then determining their roots numerically.

Short-term responses
To predict stomatal responses on hourly and daily timescales, we follow a two-step procedure. First, we find the acclimated photosynthetic capacities using the multivariate optimization described above, driven by a 7 d rolling mean of the soil water potential. Once the acclimated J max and V cmax are known, A, g s and χ can all be expressed in terms of Δψ alone. We again use the net profit in equation (1) to optimize Δψ. In this case, we determine A as the minimum of the carboxylation-limited rate A c and the light-limited rate A j . Also, since J max is fixed, the term αJ max becomes constant and can thus be ignored during optimization.

Interpretation of costs
The photosynthetic costs consist of the costs incurred by maintaining photosynthetic capacities, including the regeneration of RuBP.
Since the two photosynthetic capacities are coordinated, these costs are assumed to be proportional to J max . The hydraulic costs include (1) the construction and respiration costs of stem and leaf tissues, (2) the costs of maintaining osmotic potential and (3) the prospective costs of hydraulic failure. Since these costs are difficult to quantify through mechanistic arguments, we have taken a phenomenological approach and used the expression Δψ 2 after assessing various alternative cost expressions ( Supplementary Fig. 6). A cost expression that is quadratic in Δψ has also been adopted previously 23 . Sensitivity of our model predictions to the two cost parameters α and γ is shown in Supplementary Fig. 5.

Testing the model with data
We use published data from experiments conducted with 18 species, in which plants were grown in greenhouses under controlled conditions and subjected to progressive soil drought; values of A and g s (and sometimes also of Δψ) are reported for different values of predawn leaf water potentials, which are indicative of the soil water potential in the plant's rooting zone. The dataset has been previously assembled using tables and digitized figures from published literature as detailed in ref. 43 ,  71 For each species, data on gas exchange for different values of predawn leaf water potential were obtained from ref. 43 that originally compiled them from the sources listed in this table. Three hydraulic traits and two cost parameters were estimated using these data. A value of ∞ for drought duration means that soil water potential was maintained at each value for a long time during the experiment, allowing sufficient time for acclimation. To reduce the degrees of freedom in parameter estimation, we set b = 1 for all species, except for Helianthus annuus for which b was estimated to be 1.4.
Article https://doi.org/10.1038/s41477-022-01244-5 which we expanded to include Δψ measurements. Table 1 lists the sources and the data can be found in Supplementary Datasets 1 and 2.
In some experiments, each value of soil water potential was maintained for a long duration, so that photosynthetic capacity could acclimate (species with drought duration =∞ in Table 1). For such species, we use the multivariate optimization model as described above (equation 1).
In other experiments, the progression of drought occurred at a natural rate, ranging 12-60 d (Table 1). For such species, we use the two-step procedure outlined above to obtain the instantaneous values of the assimilation rate and stomatal conductance.

Results
We show that across 18 species, our model correctly predicts photosynthetic responses to the environment. We also show that model-predicted hydraulic strategies for the species in our dataset are consistent with widely observed empirical patterns.

Photosynthetic responses to soil moisture
Our model correctly predicts the variation in assimilation rate (A), stomatal conductance (g s ), leaf-internal-to-external CO 2 ratio (c i :c a , or χ) and soil-to-leaf water-potential difference (Δψ) in response to soil-moisture availability (ψ s ; Fig. 2). Specifically, the shapes of these dependencies closely resemble those observed during experimental drought: Fig. 3 shows predicted and observed responses for two Eucalyptus species from contrasting habitats, and Supplementary  Fig. 1 shows the corresponding responses for all 18 species. Moreover, cross-validation analysis shows that our model generalizes to out-of-sample soil-moisture conditions (Supplementary Table 1). Empirical studies report that photosynthetic capacity (V cmx and J max ) declines in response to developing soil drought 44,45 . A unique feature of our model is its ability to predict these responses correctly, qualitatively in accordance with these studies (Fig. 3e,f). Since χ depends on both J max and g s , correct predictions of χ require predicting both quantities correctly. Therefore, a close match between predicted and observed values of χ (Fig. 2c) provides further quantitative validation of the photosynthetic capacity predicted by our model.

Photosynthetic responses to vapour pressure deficit
Our work builds on the principles introduced by ref. 28 , and thus inherits the capacity to accurately predict 32 photosynthetic responses to temperature, atmospheric CO 2 and light intensity (Supplementary Fig. 2c-e). Furthermore, by explicitly accounting for plant   Fig. 2a,b). The functional relationship between logit(χ) and log(D) predicted by our model shows a close match with observations (Box 1, point 3).
In particular, our model predicts this relationship to be linear, with a median slope value of −0.697 and 5%-95% quantile range of (−0.75, −0.67) (Fig. 4a). These predicted slope values are well within the confidence interval reported in the literature 46 . Also, we find that this slope is negatively correlated with ψ 50 (Fig. 4b), such that species with highly negative ψ 50 have less negative slope values. Since earlier datasets were dominated by temperate evergreen species, this could explain why a slope value of −0.5 predicted by previous models was supported by such datasets. We offer our predicted correlation between the slope and ψ 50 as an empirically testable prediction for future studies.

Predicted hydraulic strategies match observations
In this section, we compare several widely observed empirical patterns among plant hydraulic traits with the corresponding model-predicted patterns. This qualitative comparison allows us to validate our model at an even deeper level.
First, we compare the distribution of the model-predicted degree of anisohydricity for the 18 analysed species (Box 1) to an empirically observed distribution obtained from a recently compiled database on 102 species across the globe 47 . For each species, the degree of anisohydricity is determined by the slope of the relationship between the water potential in the leaf (ψ l ) and in the soil (ψ s ), measured at low ψ s (slope <1 for isohydric, =1 for isohydrodynamic and >1 for anisohydric species). The observed global distribution of these slopes peaks at approximately 1, suggesting that the global majority of species follow the isohydrodynamic strategy 47 . The corresponding distribution  Supplementary Fig. 1. a-f, Predicted responses (lines) and observed responses (points) to decreasing soil water potential (ψ s , measured as predawn leaf water potential): assimilation rate A (a), stomatal conductance g s (b), leaf-internal-to-external CO 2 ratio χ (c), soil-to-leaf water-potential difference Δψ (d), carboxylation capacity V cmax (e) and electron-transport capacity J max (f). Eucalyptus pilularis (blue lines and squares) typically occupies warm and humid coastal areas in eastern Australia, whereas Eucalyptus populnea (green lines and triangles) typically occupies semi-arid interior regions of eastern Australia. Since both species were grown in the same greenhouse during the experiment, their contrasting responses reveal genetic adaptations to their native environments. For both species, progressive drought was experimentally induced over 12 d, resulting in a fast instantaneous response of stomatal conductance in combination with a slow acclimating response of photosynthetic capacity. Our model predictions readily account for both responses.
Article https://doi.org/10.1038/s41477-022-01244-5 predicted by our model lies within the observed distribution (Fig. 5a). Similarly, the predicted distribution of typical operating water potentials (ψ l at ψ s = 0) also closely matches the corresponding empirically observed distribution (Supplementary Fig. 3). Second, we compare the relationship between model-predicted plant hydraulic vulnerability and empirically observed turgor loss point for a subset comprising 7 of the 18 analysed species for which these empirical data were available. Empirical observations show that the turgor loss point lies between the point where the leaf loses 50% conductivity (ψ 50 ) and the point of stomatal closure (ψ g88 ) 48 . Our model predictions are consistent with this observation for most of those 7 species (Fig. 5b).
Third, we compare how model-estimated plant hydraulic vulnerability (ψ 50 ), a measure of hydraulic safety, covaries with model-estimated plant hydraulic conductance (K p ), a measure of hydraulic efficiency. Global data reveal a trade-off between safety and efficiency, that is, no plants score high on both, but only a weak correlation between them, that is, many plants score low on both. Consistent with these observations, we find only a weak correlation between model-estimated values of ψ 50 and K p , with a few species having low values of both traits, but no species having high values of both (Fig. 5c).
Fourth, we compare how the model-estimated stomatal closure point (ψ g88 ) relates to empirically observed xylem hydraulic vulnerability (ψ 50x ) for the 18 analysed species. Empirical observations show that stomatal closure occurs before the onset of substantial xylem embolism 37,49-51 , which is probably an adaptation to prevent plant mortality during drought 4 . At the same time, the minimum water potential experienced by the leaves (ψ min ) is close to the water potential at which the xylem loses 50% conductivity (ψ 50x ), leading to extremely low hydraulic safety margins 52 . Both of these observations are matched by our model, confirmed by a close correlation between ψ g88 (which is a proxy of ψ min ) and ψ 50x (Fig. 5d).

Discussion
We have presented an analytical trait-based optimality model, unifying plant photosynthesis and hydraulics to predict the stomatal responses and biochemical acclimation of plants to changing hydroclimates. Consistent with widely observed empirical patterns and benchmarked with experimental data available for 18 different plant species, our model correctly predicts the stomatal and photosynthetic responses to soil drought and the dependencies of photosynthesis on vapour pressure deficit, temperature, light intensity and CO 2 .

Comparison with other stomatal optimization models
To our knowledge, our model is probably the first to combine semi-analytical simplicity with physiological realism to predict the simultaneous stomatal and biochemical responses of plants to the environment, including water stress. Our approach has four key strengths: (1) the multivariate optimization used in our model allows for predicting the observed s-shaped decline of V cmax in response to drying soil from first principles, in contrast to most other models that require specification of V cmax as a species-specific trait; (2) explicit inclusion of plant hydraulics in our model enables separating the stomatal responses to atmospheric drought and soil-moisture availability, which could be used for improving remote-sensing estimates of GPP; (3) being based on optimality principles and specified with just two latent parameters, our model can be expected to perform well even under out-of-sample environmental conditions, that is, those outside the domain of environmental factors used for calibration, such as elevated CO 2 levels; and (4) enabling a semi-analytical solution, our model is computationally efficient and can thus be readily incorporated into existing vegetation modelling frameworks. A quantitative model intercomparison between stomatal models including ours would provide valuable insights into the relative accuracy and strengths of different stomatal optimization frameworks and would thus be an interesting direction for future research.
Leaf photosynthesis is known to be jointly constrained by stomatal and non-stomatal limitations. A vast majority of photosynthesis models account only for stomatal limitations, where stomatal conductance is optimized to maximize photosynthetic gain. Non-stomatal limitations, such as the constraints imposed by leaf mesophyll, photosynthetic capacity and sugar transport, have received attention only in the most recent stomatal models. Even such models account for them using a pre-determined functional response, in which mesophyll conductance 53 or photosynthetic capacity 24,53,54 is scaled in a prescribed way with stomatal conductance. By contrast, our model can account for non-stomatal limitations without the need to specify photosynthetic capacity a priori. This is especially important when applying it to out-of-sample environmental conditions and to species for which empirical estimates of photosynthetic capacity are not available.
Almost all current models of stomatal optimization focus on water transport through the xylem. By contrast and in line with growing evidence 37,40-42 , we hypothesized that the outside-xylem segments of the hydraulic pathway (in the leaves and roots) together form the hydraulic bottleneck of the plant. The hydraulic-segmentation hypothesis 55 states that expendable organs such as leaves and fine roots act as a 'safety valve' by losing conductivity and driving stomatal closure before the onset of fatal xylem cavitation. Therefore, we hypothesized that our model-based estimates of ψ 50 would arise from leaves and roots and thus be less negative than the corresponding empirically observed values for xylem (ψ 50x ). Our findings are consistent with these hypotheses. We offer the fitted values of these traits (ψ 50 and K p ) as testable predictions that can be validated by explicitly measuring these traits for the leaves and roots of the corresponding species.

Model assumptions and extensions
Under extreme hydroclimatic conditions, such as extremely dry or flooded soils, or extremely low atmospheric CO 2 levels, our model predictions may deviate slightly from observations. For example, most species in our dataset show an increase in the leaf-internal-to-external CO 2 ratio (χ) after stomatal closure (Fig. 3). The model-predicted χ does  46 (their reported mean and confidence interval is shown by the green line and green region, respectively). It is significantly different from −0.5 (orange line; a one-sample t-test shows a predicted mean slope of −0.7 and a 95% confidence interval of (−0.72, −0.68)). For each species, we calculate the predicted slope by varying vapour pressure deficit in the range 5-5,000 Pa while keeping other environmental parameters constant (at values reported in the respective experiments, with ψ s = 0) and using fitted trait values (Table 1). b, This slope is correlated with the ψ 50 (black points are species-specific values and the blue line is a linear regression line), with more negative slopes observed for species with less negative ψ 50 (drought avoiders). This could be a reason why earlier datasets supported a slope value of −0.5, as such datasets were often dominated by temperate evergreen species, which are typically characterized by highly negative values of ψ 50 .
Article https://doi.org/10.1038/s41477-022-01244-5 not increase, asymptotically approaching a constant value instead ( Fig. 3 and Supplementary equation 18). The increase in χ is due to a build-up of CO 2 in the leaf, which can happen via two mechanisms: (1) if dark respiration continues even after assimilation has ceased, or (2) if assimilation and respiration decline together due to reduced photosynthetic activity while CO 2 continues to 'leak in' through the leaf cuticle 56 . Future research could identify which of these processes is observed in leaves: since the source of CO 2 is different in the two mechanisms (plant metabolism or ambient air), they can be distinguished on the basis of whether the build-up is detected also in δ 13 C measurements. Our model could be extended to include these mechanisms, but this would not affect the predicted assimilation rates because these mechanisms are relevant only after stomatal closure.
For simplicity, we have assumed a single effective hydraulic pathway characterized by effective traits summarily describing water flow through the different segments (leaves, xylem and roots) of the pathway. However, a more realistic extension of our model could readily be developed by explicitly modelling water transport through each segment. In Supplementary Information section 2, we present a derivation of an extended model accounting for a multi-segment pathway.
Such an extension would be particularly useful for resolving the roles of roots and leaves in stomatal control and drought survival.
To avoid making the model too complex and parameter-heavy and to enable a semi-analytical solution, we have neglected two physiological details, inclusion of which offers promising directions for further research: (1) the leaf's energy balance and (2) hysteresis in the conductivity-soil-moisture relationship. First, when soil dries, reduced transpiration raises the leaf temperature, which in turn affects the temperature-dependent photosynthetic parameters and the dark-respiration rate. Second, the recovery of root and stem conductivity after rehydration is typically slower than the corresponding loss of conductivity during dehydration 42 . This leads to hysteresis in the relationship between conductivity and soil water potential. The speed of recovery is especially hampered under extreme drought when xylem becomes cavitated because recovery from cavitation requires embolism refilling or growth of new xylem. In our model, the relationship between soil water potential and conductivity is described by a vulnerability curve P(ψ) that does not include hysteresis. To account for hysteresis, one may use a separate vulnerability curve during recovery or use a hysteretic submodel for conductivity 39 .  47 ). b, Consistent with empirical observations, the observed turgor loss point (thick green line) lies between the model-predicted water potential at 50% loss of plant conductivity (ψ 50 ; black line) and the model-predicted water potential at 88% stomatal closure (ψ g88 ; brown line). c, Plant hydraulic conductance (K p ) is weakly negatively correlated with ψ 50 , with no species having high values of both traits, implying a weak safety-efficiency trade-off in line with empirical observations. d, When leaf water potential is at ψ g88 , the observed loss of xylem conductivity is typically less than 50% (implied by observed xylem hydraulic vulnerability ψ 50x being less than model-predicted ψ g88 ), which means that plants close their stomata before the onset of substantial xylem embolism. Furthermore, the difference between the regression line (black) and the 1:1 line (grey) is low, implying that the hydraulic safety margin ψ 50x − ψ g88 is small on average. Sources for values of ψ 50x are given in Supplementary Table 2. Closed circles indicate species for which γ was estimated using data on Δψ, whereas open circles refer to species for which such data were not available and for which we therefore used an average value of γ estimated for the respective plant types.
Article https://doi.org/10.1038/s41477-022-01244-5 Further research focusing on drought and rewetting experiments can help generate data for robustly parameterizing and validating such a model. Plants respond to increasing water stress on multiple timescales and at multiple scales of organization (leaves, whole plants and even stands). Our model accounts for leaf-level responses on daily and weekly timescales, capturing the role of stomatal closure in preventing damagingly negative stem water potentials. Scaling responses from the leaf level to the whole-plant level requires considering the distribution of leaves within the plant canopy as well as the total canopy area: this could be achieved by embedding our model within a model of plant canopies. Such a model could then be extended further to yield optimality-based first-principles models for the plant level 57 , accounting for responses beyond the point of stomatal closure, such as the shedding of leaves on a monthly timescale, which occurs either in full to prevent loss of water through cuticular tissue 4 or in part to reduce transpiration demand and continue photosynthesis 58 . Similarly, changes to the characteristics and architecture of the transpiration pathway, as reflected in stem traits 59 , could be modelled on yearly to decadal timescales. Trait adaptation on centennial or millennial timescales could be modelled by embedding our leaf-level optimality theory into evolutionary models [60][61][62] .
In our model, on the shortest timescale (minutes to days), plants may optimize leaf water potential for a fixed (acclimated) photosynthetic capacity. On the timescale of weeks, plants additionally adjust their photosynthetic capacities. In principle, the weekly timescale can be modelled either with nested optimization (that is, optimizing daily, or sub-daily, stomatal conductance for a given photosynthetic capacity, and optimizing weekly photosynthetic capacity by maximizing the total profit over a week), or with simultaneous optimization (that is, optimizing both variables together by assuming a constant environment during the week, representative for mean daytime conditions). In this work, we have taken the latter approach for theoretical and computational simplicity, but the alternative nested approach is worth exploring in future work.

Implications for global vegetation modelling
Taking advantage of the increasing availability of global data on plant traits, our model can be applied at the global scale by making a few additional assumptions. For the species in our dataset, we find that the photosynthetic unit cost α lies within a narrow range of values (0.08-0.12), which could therefore be treated as a constant. Furthermore, variations in the hydraulic unit cost γ are primarily driven by differences between plant types, with relatively lesser variability within plant types, suggesting that γ could also be treated as a constant within plant types (Supplementary Fig. 4). To infer the global distribution of α and γ, our model can be used in a Bayesian framework on global data on gas-exchange measurements. When plant-level hydraulic traits ψ 50 and K p are not available, these could be derived from other widely measured traits on the basis of observed patterns of functional coordination among plant organs. As a starting point, we have shown that K p is weakly correlated with ψ 50 (Fig. 5c) and specific leaf area (Supplementary Fig.  4b). Further studies could test for similar correlations with leaf vein density and root mass fraction, which are respectively expected to affect leaf and root conductivities.
Accurate models of plant photosynthesis are crucial for improving projections of the global carbon and water cycles, because photosynthesis and transpiration by terrestrial plants account for 56% and 30% of the global fluxes of carbon dioxide and water, respectively 63,64 . It is especially important to develop models that can generalize to new climatic conditions, because the projected increase in the frequency and intensity of droughts may lead to unprecedented climatic conditions in future. The inclusion of plant hydraulics into vegetation models has been shown to improve predictions of global productivity and evapotranspiration [17][18][19][20][21] , as well as predictions of the spatiotemporal diversity of vegetation 65 . Spearheading the development initiated by these studies, our model is ideally suited for being embedded into global vegetation models: by accounting for biochemical acclimation, plant hydraulics and photosynthetic trade-offs through optimality principles, our model can extend to new species and new environmental conditions with a raised degree of confidence. Furthermore, accounting for photosynthetic and hydraulic costs is bound to yield more accurate estimates of the energy spent on resource acquisition and, consequently, of the resources available for growth and reproduction. Therefore, embedding our model of photosynthesis and hydraulics into a demographic model of plant communities can help improve the scaling of assimilation and transpiration from the leaf level to the whole-plant level, and even from plants to communities, thus paving the way for more accurate and robust land-surface models.

Methods
Our model consists of three key components or modules, corresponding to the three principles and hypotheses: a water-transport module to account for plant hydraulics and water balance, a photosynthesis module to account for photosynthesis and the photosynthetic-coordination hypothesis, and a profit-maximization module to implement the optimization. Here we describe the equations used for each module, as well as our strategy for model calibration. Full derivations of the equations can be found in Supplementary Information section 1.

Water-transport module
We model water transport using Darcy's law applied to small cross-sections of the hydraulic pathway ( Supplementary Information  section 1.1.4). In principle, our model of water transport can explicitly represent multiple segments ( Supplementary Information section 2), but for simplicity, we represent the entire pathway by a single 'effective segment' with traits K p , ψ 50 and b. Thus, our hydraulic model is mathematically similar to the model for xylem water transport described in ref. 24 , but the effective hydraulic traits in our model correspond not necessarily to the xylem but to the whole plant. As the outside-xylem segments (in leaves and roots) are the hydraulic bottlenecks of the plant for many species 37,[40][41][42] , the traits modelled here probably correspond to these segments.
The conductivity κ of any cross-section of the pathway declines as the water potential becomes increasingly negative. This decline in conductivity is phenomenologically described by a so-called vulnerability curve P(ψ), such that κ(ψ) = κ(0)P(ψ). The vulnerability curve is typically described by a Weibull function with two parameters: the water potential ψ 50 at which 50% conductivity is lost and a shape parameter b that determines the sensitivity of conductivity loss to water potential, (2) Water potential drops continuously along the hydraulic pathway, from ψ s in the soil to ψ l at the leaf surface, with a continuous decline in conductivity along the pathway. The volumetric flow of water per unit leaf area in the pathway, Q, is therefore described by a differential equation (see Supplementary Information section 1.1.4 for derivation), which can be solved for Q as follows, where κ p is the effective conductivity of the whole plant per unit leaf area, L is the effective length of the hydraulic pathway and η is the dynamic viscosity of water. The composite term κ p /L = K p is the whole-plant conductance per unit leaf area (Supplementary Information section 1.1). To keep the number of parameters low to prevent overfitting of the model to data, we use equation (3) to model water flow. In Supplementary Information section 2, we propose an extension Article https://doi.org/10.1038/s41477-022-01244-5 of this model, deriving an expression for Q by explicitly accounting for segments of the hydraulic pathway in the roots, stem and leaves. The water-balance principle states that the atmospheric demand for water imposed by vapour pressure deficit at the leaf surface matches the supply of water from the soil via the stem and leaf segments of the hydraulic pathway. The transpiration rate at which water vapour diffuses out of the leaf into the atmosphere is given by (see Supplementary Information section 1.1.6 for derivation) where g s is the stomatal conductance and D is the atmospheric vapour pressure deficit divided by the atmospheric pressure. This rate E equals the rate Q at which water enters the leaf according to equation (3), which allows us to calculate g s by solving

Photosynthesis module
The assimilation rate A is calculated from the Farquhar-von Caemmerer-Berry biochemical model 34 31 .
The photosynthetic-coordination hypothesis states that under typical daytime conditions, assimilation operates at the point of co-limitation, such that the carboxylation-limited and light-limited assimilation rates are equal. With this assumption, the co-limited assimilation rate can be written as ( Supplementary Information section 1.2 where J is the effective electron-transport capacity, which increases with light availability I abs but saturates due to limitation by the leaf's intrinsic maximum electron-transport capacity J max , Here, c a is the atmospheric CO 2 concentration, χ is the ratio of the leaf-internal and external CO 2 concentrations (c i :c a ), Γ * is the light-compensation point, K M is the Michaelis-Menten coefficient for C3 photosynthesis, ϕ 0 is the quantum yield efficiency, I abs is the absorbed photosynthetically active radiation and b r is the ratio of dark respiration to carboxylation capacity (dark respiration is assumed to be proportional to carboxylation capacity, that is, R d = b r V cmax ). Temperature dependencies of Γ * and K M are modelled according to ref. 31 . The ratio b r also has a weak dependence on temperature 66 , which we have ignored in this work. Variation in J max in response to light and water availability (by optimization) implies a coordinated variation in both carboxylation and electron-transport capacities.

Profit-maximization module
We assume that plants maximize net assimilation (or profit, F) defined in equation (1). Without loss of generality, we also assume that the unit benefit of assimilation is one, that is, α and γ represent the ratios of the unit costs to unit benefits of assimilation. To optimize equation (1), we express all quantities in terms of the two independent variables χ and Δψ and set the gradient of the profit function to 0. This can be done analytically (Supplementary equation 16). However, except in the special case of strong J max limitation, the roots of the gradient must be found semi-analytically ( Supplementary Information section 1.3.2). Solving for optimal χ * and Δψ * in turn allows us to predict the optimal photosynthetic capacities (V * cmax and J * max ), stomatal conductance ( g * s ) and CO 2 assimilation rate (A * ).

Parameter estimation and model testing
We drive the model with environmental variables (temperature, vapour pressure deficit, light intensity and CO 2 ) as specified in the experimental studies. Other parameters used in the model are as follows: ϕ 0 = 0.087, b r = 0.002. In the case of instantaneous responses, we use the daily maximum light intensity under growth conditions to calculate the acclimated response, and saturating light intensity (as specified in the studies) to calculate the instantaneous response, so as to mimic the conditions present during LiCor measurements.
Since measurements of effective whole-plant hydraulic traits are not readily available, we treat them as model parameters and estimate them along with the two cost parameters. Values of Δψ were reported for 6 of the 18 species from the same drought experiments. We supplement Δψ observations with typical values reported in the literature 67 for two additional species (Pseudotsuga menziesii and Olea europea var. Meski). For such species (for which measurements of Δψ are available), we calibrate five parameters (α, γ, ψ 50 , b, K p ) by minimizing the sum of squared errors (E r ) between predicted and observed values of A, g s , χ and Δψ, defined as where i represents different values of ψ s , E[] denotes the mean value, and variables with tilde (for example, χ ) represent observations. From this calibration, we obtain the mean estimated value of γ for each plant type. For all other species (for which measurements of Δψ are not available), we use this mean value of γ and estimate the remaining four parameters. To further reduce the degrees of freedom in model parameterization, we fix the value of b = 1 for all species, except for H. annuus for which b also had to be estimated. For each species for which data on Δψ are available, we evaluate model performance using fivefold cross-validation (or leave-one-out cross-validation when data points are limited).

P-hydro R and C++ packages
R code to run our model (P-hydro) is provided as an extension of the 'rpmodel' package (the version used in this paper is archived at https:// github.com/jaideep777/rpmodel/releases/tag/v1.0.0h), with options to use the acclimating and instantaneous responses. A C++/Rcpp version is also provided for potential integration with vegetation models (https://github.com/jaideep777/phydro). These packages provide the option to either use the semi-analytical solution derived in this work, or to directly optimize the profit function numerically. The numerical implementation also allows for quick extension of the model with different profit and cost functions. The C++ version also allows the use of an approximate calculation of g s , which substantially improves computational speed with only a minor loss of accuracy.

Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
All data used in this manuscript are compiled from the literature. We have provided citations to publications and databases at appropriate Article https://doi.org/10.1038/s41477-022-01244-5 locations in the manuscript. The compiled database can be found in the Supplementary Information.

Code availability
R code to run our model (P-hydro) is provided as an extension of the 'rpmodel' package (the version used in this paper is archived at https:// github.com/jaideep777/rpmodel/releases/tag/v1.0.0h), with options to use the acclimating and instantaneous responses. A C++ version is also provided for potential integration with existing vegetation models (https://github.com/jaideep777/phydro). The vignettes folder of this package also contains R code to reproduce key results.
Corresponding author(s): Jaideep Joshi Last updated by author(s): Jul 8, 2022 Reporting Summary Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Portfolio policies, see our Editorial Policies and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted

Software and code
Policy information about availability of computer code Data collection No new data collection was done in this manuscript. Online figure digitizer (https://automeris.io/WebPlotDigitizer/) was used to extract data from published sources.

Data analysis
Custom code used for the analysis is available on GitHub: https://github.com/jaideep777/phydro For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Portfolio guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A description of any restrictions on data availability -For clinical datasets or third party data, please ensure that the statement adheres to our policy The following data-availability statement is included in the manuscript: All data used in this manuscript are compiled from the literature. We have provided citations to publications and databases at appropriate locations in the manuscript. The compiled database can be found in the supplementary information.