Federico Vazquez
email: federico@ifisc.uib-csic.es

Cristóbal López
email: clopez@ifisc.uib-csic.es

Justin M Calabrese
email: justin.calabrese@ufz.de

Miguel A Muñoz
email: mamunoz@onsager.ugr.es

Dynamical phase coexistence: a simple solution to the "savanna problem"

Keywords: savannas, dynamical coexistence, contact process, statistical physics, stochastic processes

We introduce the concept of dynamical phase coexistence to provide a simple solution for a longstanding problem in theoretical ecology, the so-called "savanna problem". The challenge is to understand why in savanna ecosystems trees and grasses coexist in a robust way with large spatiotemporal variability. We propose a simple model, a variant of the Contact Process (CP), which includes two key extra features: varying external (environmental/rainfall) conditions and tree age. The system fluctuates locally between a woodland and a grassland phase, corresponding to the active and absorbing phases of the underlying pure contact process. This leads to a highly variable stable phase characterized by patches of the woodland and grassland phases coexisting dynamically. We show that the mean time to tree extinction under this model increases as a power-law of system size and can be of the order of 10,000,000 years in even moderately sized savannas. Finally, we demonstrate that while local interactions among trees may influence tree spatial distribution and the order of the transition between woodland and grassland phases, they do not affect dynamical coexistence. We expect dynamical coexistence to be relevant in other contexts in physics, biology or the social sciences.

I. INTRODUCTION

Savannas are open systems that feature a continuous grass layer and a discontinuous tree layer. They appear across a wide range of climatic and ecological conditions, and are characterized by the stable, though variable, coexistence of two distinct types of vegetation, trees and grasses [Sarmiento 1984]. This coexistence is dynamic in the sense that the density of trees varies widely both in space and time, as recently confirmed by observations of cyclic transitions between empty and dense tree distributions [START_REF] Moustakas | Spacing patterns of an Acacia tree in the Kalahari over a 61-year period: how clumped becomes regular and vice versa[END_REF], Wiegand et al. 2006]. Savannas have been studied from both experimental and theoretical points of view, and have become an important subject of study in ecology. However, the origin, nature, and dynamics of savannas are not yet well understood. How is long-term coexistence of trees and grasses possible without the superior competitor taking over, as happens in other ecosystems (grasslands or woodlands)? This is one of the long standing puzzles in theoretical ecology, commonly referred as the savanna problem. A tentative answer to this question is provided by niche models in which, assuming (soil, rain, etc) heterogeneities, each life form occupies the regions for which it is a superior competitor [Walter 1971, Walker et al. 1981, Walker and Noy-Meir 1982].

This type of solution is conceptually unsatisfactory and is, anyhow, not supported by recent empirical observations [START_REF] Jeltsch | Ecological buffering mechanisms in savannas: A unifying theory of long-term treegrass coexistence[END_REF], Higgins et al. 2000[START_REF] Bond | What limits trees in C4 grasslands and savannas? Annual Review of Ecology[END_REF], Scholes and Archer 1997]. Demographic (stochastic) bottleneck models [START_REF] Sankaran | Tree-grass coexistence in savannas revisited -insights from an examination of assumptions and mechanisms invoked in existing models[END_REF]] invoke stochastic explanations that rely on demographic and environmental fluctuations to generate dynamical heterogeneities [Chesson 2000, D'Odorico et al. 2007, Jeltsch et al. 2000, Jeltsch et al. 1996, Higgins et al. 2000, van Wijk and Rodriguez-Iturbe 2002, Meyer et al. 2007].

For example, the storage effect hypothesizes that long-lived organisms promote species coexistence in communities with birth rate variability, so that in an environment which is frequently adverse, a long life span buffers trees against extinction [Chesson 2000]. Other, similar in spirit, buffering mechanisms have also been proposed [START_REF] Jeltsch | Ecological buffering mechanisms in savannas: A unifying theory of long-term treegrass coexistence[END_REF]]. Alternative deterministic A c c e p t e d m a n u s c r i p t 2 descriptions, which do not rely on the assumption of rooting niche separation, have also been considered in [Scheiter 2007].

Our goal is to build a minimalistic savanna model that allows us to assess the contributions of the following features to long-term tree-grass coexistence: i) Variable weather conditions [START_REF] Jeltsch | Tree spacing and coexistence in semiarid savannas[END_REF], Higgins et al. 2000, van Wijk and Rodriguez-Iturbe 2002],

ii) Mean annual precipitation, which has been reported to enhance and limit the maximum tree-cover [START_REF] Sankaran | Determinants of woody cover in African savannas[END_REF], Sankaran et al. 2008, Bucini and Hanan 2007], iii) Adult tree longevity [START_REF] Higgins | Fire, resprouting and variance: a recipe for grass-tree coexistence in savanna[END_REF],

and iv) Positive and negative local density-dependent tree interactions [START_REF] Jeltsch | Tree spacing and coexistence in semiarid savannas[END_REF], Meyer et al. 2007, see Calabrese et al. In press].

Other ingredients to be studied are fire and dispersal limitation. Some existing demographic savanna models, which include these elements along with many others, can reproduce the main traits of real savannas [START_REF] Jeltsch | Ecological buffering mechanisms in savannas: A unifying theory of long-term treegrass coexistence[END_REF], Jeltsch et al. 1996, Higgins et al. 2000, Sankaran et al. 2004, van Wijk and Rodriguez-Iturbe 2002, Meyer et al. 2007]. However, these models do not clarify which ingredients are necessary to produce long-term coexistence and which are superfluous.

The gist of our minimal modeling is as follows. We start from a cellular automata defined on a square grid, in which each site can be occupied either by one tree or grass, and that follows the dynamical rules of the standard contact process (occupied sites are trees and empty ones are grasses; see next section). As in the standard contact process, depending on parameter values, two phases characterize the system: grassland (grass only) and woodland (tree dominated). "Dynamical phase coexistence" appears when the system is driven by a random external driver (fluctuations in rainfall or other environmental effects) which forces the system to visit the two phases randomly in time, alternating between grassland prone conditions and woodland prone condition. We shall show that the range of dynamical coexistence is much enhanced if the age of the trees is included, so that trees can typically endure harsh conditions. In any case, coexistence appears in a broad region of parameter space without fine-tuning of parameters. This dynamical coexistence is not indefinite as the environmental fluctuations will eventually lead the system to the grassland or absorbing state, which is characterized by the complete absence of trees.

Unlike previous work on environmental fluctuations and tree-grass coexistence

[van Wijk and Rodriguez-Iturbe 2002, [START_REF] Jeltsch | Tree spacing and coexistence in semiarid savannas[END_REF], Higgins et al. 2000], we focus here on characterizing the timescales over which dynamical coexistence is robust. Our approach consists of studying stability by means of the mean life-time, that is, the mean time that takes the savanna to reach the final absorbing state. We show that this time diverges with system size N as N α , and can be enormous for even moderately sized systems [Leigh 1981].

Finally, let us remark that the generic phase coexistence observed in our model, i.e. coexistence occurring in a broad region of parameter space, is rarely found in other non-equilibrium model/systems and, in most cases, the domain of coexistence is typically small (see [Muñoz et al. 2005] and references therein). Additionally, given that the phenomenon of coexistence also happens in biological, physical and social sciences, the mechanism proposed in this paper is expected to be useful in many other contexts.

II. MODEL

Consider the simple contact process (CP) [Harris 1974, Hinrichsen 2000]. Each node (i, j) of a two-dimensional square lattice can be either occupied z i,j = 1 (tree) or vacant z i,j = 0 (grass). The dynamics is as follows: a tree is randomly selected, and it is removed from the system with probability d, otherwise, with probability b, it generates an offspring, which is placed at a randomly chosen nearest neighbor (n.n) provided it was empty (i.e. short-range seed-dispersal).

Similar base assumptions underpin a number of CA savanna models including those of Menaut et al. (1990), Hochberg et al. (1994), [van Wijk and Rodriguez-Iturbe 2002], Kefi et al. 2007, and [see Calabrese et al. In press].

Even much more elaborate models such as those by [START_REF] Jeltsch | Tree spacing and coexistence in semiarid savannas[END_REF]] and [START_REF] Higgins | Fire, resprouting and variance: a recipe for grass-tree coexistence in savanna[END_REF]] have a similar core.

Thus the contact process (together with its extensions) is a widely used, though rarely explicitly acknowledged (but see [see Calabrese et al. In press]), basis for a savanna model. Here, we take full advantage of the body of results that has developed around the contact process, and in doing so we place our savanna model on firm theoretical ground. We then to build on this base in a systematic way and assess the contributions of the additional model components detailed below to tree-grass coexistence. Our goal is to identify a minimal set of processes that promote robust, longterm tree-grass coexistence, and we make no attempt to build a model that is a quantitatively accurate representation of any particular savanna. We therefore intentionally leave out many details considered in site-specific savanna models (e.g. [START_REF] Jeltsch | Tree spacing and coexistence in semiarid savannas[END_REF], Higgins et al. 2000, Meyer et al. 2007]).

As is usually done in sequentially updated Monte-Carlo simulations, every time a tree is selected, time t is increased by 1/N (t), where N (t) is the total number of trees in the system, so that when all trees have The underlying pure model has two homogeneous phases: active (woodland) and absorbing (grassland), separated by a critical point (the transition could also be discontinuous as in Fig. 4). The control parameter fluctuates in time and shifts from one phase to another (droughts and rainy periods). This, combined with long-living trees, prevents the system from reaching a homogeneous steady state.

been selected once on average, t is increased in one unit, corresponding to one Monte Carlo (MC) step or one year. The dependence of the timestep on N (t) simply means that the more trees in the system the more individual updatings need to be done per each MC step. Fixing d, a phase transition appears at some critical value b c . For b > b c the system is in the active (the system has a non-vanishing density of trees and can dynamically evolve) phase -woodland -, while for b < b c it is in the absorbing -grassland -phase with zero tree-density (as schematically illustrated in Fig. 1). The reason for this last name is that if the system reaches this absorbing (i.e. only-grass) state, the situation is irreversible, it remains indefinitely trapped in it. Since there is no spontaneous generation of trees, the dynamics ceases if there are only grasses. This type of transition is commonly termed an absorbing phase transition [Hinrichsen 2000]. Now, we introduce the following two extra ingredients to model savannas:

i) Fluctuating external conditions: We assume that the birth probability depends on external conditions (mostly annual precipitation, but also fires, etc [START_REF] Sankaran | Determinants of woody cover in African savannas[END_REF]), so b becomes a timedependent random variable (see the zigzaging line in Figure 1). To account for possible temporal correlations in weather conditions we take a time-correlated (colored) noise as follows. With probability q a new value of b is extracted from an uniform random distribution in [b min , b max ] at each MC step, otherwise (with probability 1q) b is kept fixed. Here, we take b min = 0, while b max is the control parameter. This models weather cycles of typical length 1/q. In a similar way, one could include more periodic weather oscillations, as those induced by "El Niño", which lead to similar results.

ii) Age: Field studies reveal that the mortality distribution of some savanna trees is consistent across years and climatic regimes [Moustakas et al. 2006], and that they have long lifespans. To model this we define trees with an intrinsic age-variable, a(i, j) measured in years. The death probability is taken to be age-dependent: d → d(a(i, j)). In particular, a random number, η, is extracted from a Gaussian distribution of mean a m and variance σ (typically, a m = 100 and σ = 20). If a(i, j) ≥ η then the selected old tree is removed; otherwise nothing happens.

Some other effects, such as density-dependent demographic rates, can also be easily implemented in the model.

Density-dependence: Negative and positive local density-dependent death probabilities account for tree-tree competition and facilitation, respectively. Both of these effects have been reported to act in savannas [see Calabrese et al. In press]. To model competition between a young tree at site (i, j) (tree with age a(i, j) ≤ a est , where a est is the establishment age) and its neighbors, we increase its death rate as a function of the number of its nearest neighbors trees N N (i, j) = (k,l) z k,l , where the sum is restricted to nearest neighbors sites of (i, j). Thus, we take d(i, j) = 1exp [-N N (i, j)]. Contrarily, to model strong facilitation we consider d(i, j) = 1 for N N (i, j) = 0 and 1, and d(i, j) = 0 for N N (i, j) ≥ 2, that is, trees born in sites with only a few occupied neighbors die with certainty, otherwise they survive.

Fire: while fire is widely recognized as an important aspect of the ecology of many mesic to wet savannas, is clearly of secondary importance relative to the role played by moisture availability [START_REF] Sankaran | Determinants of woody cover in African savannas[END_REF], Sankaran et al. 2008, Bucini and Hanan 2007[START_REF] Bond | What limits trees in C4 grasslands and savannas? Annual Review of Ecology[END_REF]. Over the long time scales which are our focus here, the dynamics of tree cover and the coexistence of trees and grasses will be dominated by fluctuations in climatic conditions. Over shorter time scales and when the details of tree spatial pattern are the main focus, fire should certainly be considered in more detail than we have here. However, given our long-term focus and our goal of identifying the primary driver of tree-grass coexistence on these time scales, its main effect can just be accounted for by the fluctuating parameter b; years with (without) big fires will have smaller (larger) b.

Dispersal limitation The role of dispersal limitation in savannas is not completely resolved yet. However, it seems that there is some tendency for high seed densities near trees [Barot et al 1999, Witkowski and Garner 2000] and,

III. MODEL ANALYSIS AND RESULTS

We first analyze the role of each ingredient we have added to the basic contact process separately. Taking a fluctuating b and a fixed death rate, i.e. switching off the age effect, the system shifts randomly between the tendencies to be in the active (tree density larger than zero) and in the absorbing (zero tree density) phase of the underlying pure model (notice the zigzagging trajectory in Fig. 1). For b max > b c , the system hovers around its critical point, while, if b max < b c (resp. b min > b c ) the fluctuating system is in the absorbing (active) phase, i.e. it is a grassland (woodland). In principle, if the time series of b happens to be adverse (i.e. b < b c ) for a sufficiently long time interval, any finite system falls into the absorbing state; i.e. the system has variability but little resilience to long adverse periods.

The effect of weather correlations is as follows: for 0 < q << 1 the birth rate is constant for long periods and the system typically jumps, every 1/q years, from a pure CP homogeneous state to another one and, therefore, when b takes a value smaller than b c it falls ineluctably into the absorbing state. For q ≈ 1, the rate of variation of b is very fast and the system does not have the time required to relax to any pure CP steady state, and reaches instead an averaged density value (see solid curve in Fig. 2). For realistic intermediate values (e.g. q = 0.03), the system exhibits much larger oscillations (see the dashed and the dotted line in Fig. 2) which resemble those in real savannas [START_REF] Moustakas | Spacing patterns of an Acacia tree in the Kalahari over a 61-year period: how clumped becomes regular and vice versa[END_REF]].

On the other hand, taking a fixed birth rate, q = 0, i.e. switching off the weather variability the model becomes a CP with age. For this, if trees die at a fixed given maximum age, a m (i.e. with σ = 0), the density is known to exhibit damped oscillations in time of period 2a m , and to converge asymptotically to a homogeneous stationary value [Gerami 2002] (see dashed-dotted line of Fig. 2). Our model (with σ > 0) exhibits analogous, though more variable, damped oscillations, and converges either to the absorbing state or to an active homogeneous state. Separating these two regimes there is a phase transition. In this case the system is resilient but has little variability. Obviously resilience grows with the maximum age.

As we demonstrate below, the full model, which includes both age and correlated variable rainfall, exhibits variability and resilience (see the dashed and dotted curves in Fig. 2); the larger the maximum age, the larger the resilience. The system fluctuates locally between the absorbing and the active phases, The dashed-dotted curve (q = 0, i.e. model without variability) exhibits damped oscillations. The solid line (q = 1) shows small variability, while the dashed and dotted curves are for q = 0.03 (intermediate variability). Notice that the one with smaller maximum age (am = 30; dotted curve) does not survive to an adverse period, while its analogous for large age (am = 100; dashed curve) does. but it is able to preserve "islands" of the unfavored phase in a "sea" of the dominant phase, as is required for generic phase coexistence. This is the basic mechanism of dynamical phase coexistence.

The effect of density-dependence is as follows. Computer simulations show that, the underlying pure model (q = 0 and no age), either in the absence of density dependence and in the case of competition or weak facilitation, exhibits a continuous absorbing phase transition (as schematized in Fig. 1). Instead, for strong facilitation the underlying transition can be discontinuous with a broad hysteresis loop (see Fig. 3), implying that around the transition the two dynamically coexisting phases are very different.

Note that the appearance of a discontinuous transition in the underlying pure model is interpreted in biological terms as the ex- istence of catastrophic shifts that may result in sudden dramatic changes of the savanna properties [START_REF] Scheffer | Early-warning signals for critical transitions[END_REF], Scheffer et al 2003] Fig. 4 shows the stationary tree density as a function b max (q = 0.1, a m = 100, a est = 1) illustrating the existence of an active phase and a continuous phase transition in the full model with competition. The active phase of the full model is the phase of coexistence: the two phases of the pure model (grass and wood) coexist. A very similar continuous transition is obtained in the case of strong facilitation; even if the underlying pure-model transition is discontinuous, once varying conditions are switched on, the transition between the absorbing and the tree-grass coexistence phases becomes continuous.

For illustration purposes, Fig. 5 shows snapshots of such a phase for different parameters at different times. Panels (a) and (b) correspond to the case of competition, and both have the same parameter values but look quite different, illustrating the large spatio-temporal variability. Panel (c) shows that, in the case of facilitation more compact clusters are observed (as justified by the underlying discontinuous transition). These snapshots are visually very similar to pictures of real savannas [START_REF] Moustakas | Spacing patterns of an Acacia tree in the Kalahari over a 61-year period: how clumped becomes regular and vice versa[END_REF]; notice the presence of irregularly distributed tree clusters of different sizes and shapes. Panel (d) corresponds to competition, but with slightly different parameters. We can conclude that density-dependence controls the shape of emerging clusters, but it is not an essential ingredient for dynamical coexistence.

Let us now scrutinize the stability of the active phase, where there is dynamical coexistence: how long does it persist through periods of adverse rainfall conditions?

In a single realization, independent of the value of b max and provided that b min is smaller than the critical values of the pure model, the density of trees reaches the absorbing state ρ = 0, that consists of a totally empty lattice where no more trees are born, and therefore the evolution stops. In other words, trees go extinct in a time that depends on the specific realization, due to stochastic fluctuations, and the size of the system. We have computed the mean lifetime, τ , that is, the average over many realizations of the time required to reach tree extinction, for several values of b max and different system sizes N (see Fig. 6). The value b c max 0.1045 at which τ grows as τ ∼ (ln N ) 3.68 [straight dashed line in Fig. 6(a)], agrees within error-bars with the point where the stationary density ρ s goes to zero in Fig. 4, thus we take it as the active-absorbing phase transition point. In Fig. 6(b) we observe that, remarkably, in the active phase (circles) τ diverges with N as a power law N α (with possible log corrections), where the exponent α increases with b max (see inset). That is, the mean life-time increases very rapidly with system size, and therefore one expects that real (large) savannas survive for very long periods. A similar behavior for the extinction time was found in a discrete-time population model (the lottery model) in the presence of a stochastic environment in [Chesson 1982]; waiting times for extinction increase at very different rates: power law, logarithmic and square-logarithmic. The author discusses the mathematical and biological properties of the system that give rise to each of these behaviors.

We now provide a simple analytical explanation for the power-law behavior of the extinction time with N . The relaxation time T relax that the system takes to reach the absorbing state, when b takes values below b c , scales as T relx ∼ ln(N ), because the density of trees decays exponentially fast to zero in the absorbing phase. The system dies when there is an "extreme event" consisting of a dry weather period larger than T relx , that is, when b takes a value bellow b c a number n = qT relax of consecutive times, given that b takes a new value every average time interval 1/q. Then the probability of having an extreme event is P ∼ (b c /b max ) n , so that the mean time τ to observe this event by the first time is τ ∼ 1/P ∼ (b c /b max ) -q ln(N ) ∼ N α , with α ≡ -q ln(b c /b max ), showing the algebraically long surviving times. Note that the exponent α is proportional to b max , as it is observed in the inset of Fig. 6(b).

As an illustration of the mean-life time, consider a 1000 hectare savanna, which corresponds to a 632 × 632 square lattice (assuming, say, neighboring lattice sites are separated by 5 meter distance), taking b max 0.111, which is not very deep into the active phase, then our model predicts that its expected time to extinction is of the order of 10 million years. This leads to the conclusion that the coexistence of grassland and woodland phases is stable for extremely long time periods. Furthermore, we find that the exponent α increases with the maximum tree age a m (not shown), so that age strongly stabilizes the active phase and increases the savanna mean life-time time.

Our main conclusion is that, in the infinite size limit, the dynamical coexistence of savannas is truly stable, even if, for any finite size, surviving times are necessarily finite. The corresponding extinction times increase as power-laws of system-size; this is, systems with dynamical phase coexistence are less stable than homogeneous steady states (for which there is an exponential increase of surviving times upon enlarging system size), but still they are stable.

IV. EMPIRICAL EXAMPLE

To demonstrate that this simple model is a reasonable descriptor of savanna dynamics, we fit it to a paleoecological time series of tree-grass pollen ratios, R data (t), and lake depth levels, δ(t), from the Crescent Island Crater site in central Kenya [START_REF] Lamb | Vegetation response to rainfall variation and human impact in central kenya during the past 1100 years[END_REF], Vershuren et al. 2000]. The long-term nature of this time series (> 1000 years) allows us to examine a much wider range of density fluctuations than is possible with short-term, directly observed tree density data. The cost is that these paleoecological data sets rely on proxy measures (e.g. fossil pollen for vegetation, fossil midge and diatom assemblages for lake depth) and are thus necessarily much less precise than direct measures of tree density and rainfall. We digitized figures 3B and 3C in [START_REF] Lamb | Vegetation response to rainfall variation and human impact in central kenya during the past 1100 years[END_REF]] to obtain the pollen ratio and lake depth time series, respectively, using the program Engauge V 4.1 [http://digitizer.sourceforge.net/]. As these time series were both reconstructed at the same site, over the same time interval, they give us rough proxy measures for both precipitation patterns and tree density at this site going back over 1000 years (see Fig. 7).

We ran simulations of the model with a birth probability given by the expression b(t) = γ δ(t) (where γ is a constant), that is, assuming that the amount of rainfall, and therefore b(t), is proportional to the lake's depth. We included neither age nor density effects. We also considered that the amounts of pollen from trees and grasses are proportional to the number of trees ρL 2 (with L = 100) and the area covered by grass (1ρ)L 2 , respectively, thus we took R num (t) = β ρ/(1ρ) as the tree-grass pollen ratio in the simulations, where β is a constant. Then, for given values of γ and β, we averaged over 10 independent realizations of the dynamics, starting from the same initial density of trees, and calculated the sum of squared deviations between the model and the tree-grass ratio data over the 1000 years period as

1000 i=1 R i data (t) -R i num (t) 2 .
Comparing the model to the paleoecological data from the Crescent Island Crater site, we found by numerical direct search that γ min = 0.018 and β min = 0.005 are the values that minimized the sum of squared errors. As can be seen in Fig. 7, the average time evolution of the tree-grass pollen ratio (solid line) is strongly biased by the temporal function b(t) (dashed line), and it has only a qualitative agreement with the evolution of the real ratio (dashed-dotted line). We have checked that the age of trees in the model does not affect the evolution of the ratio very much, and that single realizations are typically very similar, given that fluctuations coming from the system size are much smaller than fluctuations induced by the variation of b(t). We also note in Fig. 7 that the real-data pollen ratio also roughly follows the variation of the lake's depth. These results demonstrate that this correlation between rainfall and pollen ratio, already noted in [START_REF] Lamb | Vegetation response to rainfall variation and human impact in central kenya during the past 1100 years[END_REF], Verschuren et al. 2000], is correctly described by our very simple model. This correspondence is especially encouraging given that we have manipulated only two parameters here to fit the model.

V. SUMMARY

Here we have shown that a simple extension of the contact process is capable of providing a robust and general explanation for tree-grass coexistence in savannas. Specifically, fluctuating rainfall levels force the system to oscillate between a tree-dominated woodland (active) phase and a grass-only (absorbing) phase; a behavior we call dynamical phase coexistence. Dynamical coexistence is not permanent, but, as we have shown, it is expected to last for geologically significant periods of time. In other words, the mechanism discussed here facilitates coexistence over periods of time easily long enough to span the gaps between major disturbance events, such as ice ages, that can create and destroy savannas.

Strictly speaking, the only essential ingredient for dynamical phase coexistence is the presence of externally varying conditions. Said another way, dynamical coexistence is independent of the lifespans of trees, the degree of correlations in weather conditions, and the presence and nature of density-dependent local interactions among trees. Adding age-dependent tree death to the model greatly enhances the stability and robustness of coexistence. Adding weather correlations of intermediate magnitude increases the variability in tree cover and can help the model mimic the wide variability in tree cover that has been observed among sites with similar rainfall [START_REF] Sankaran | Determinants of woody cover in African savannas[END_REF], Bucini and Hanan 2007, Sankaran et al. 2008]. Finally, adding local treetree interactions can allow the model to reproduce a range of tree spatial patterns that have been observed in real savannas, including dense thickets of trees and more widely spaced and open configurations [see Calabrese et al. In press, Moustakas et al. 2008]. Thus these additional features affect the nature of several observable features of savannas but are not the primary drivers of treegrass coexistence. This last point merits further discussion. While several stochastic savanna models have included varying weather conditions and many other factors [START_REF] Jeltsch | Tree spacing and coexistence in semiarid savannas[END_REF], Jeltsch et al. 1998, Higgins et al. 2000, Meyer et al. 2007, van Wijk and Rodriguez-Iturbe 2002], none has, to our knowledge, clearly identified the minimal conditions that facilitate tree-grass coexistence. Our strategy of starting from the well-studied contact process and extending the model in a stepwise fashion has allowed us to identify the contribution of each model component to tree-grass coexistence and to other features observed in real savannas. We can therefore unambiguously state that fluctuating external conditions alone are sufficient to facilitate long-term (though not indefinite) treegrass coexistence.

A further advantage of our minimalistic approach is that it does not require parameter fine tuning to achieve coexistence, as has been the case in other stochastic savanna models [START_REF] Jeltsch | Tree spacing and coexistence in semiarid savannas[END_REF], Higgins et al. 2000].

Finally, unlike previous studies, we have been able to clearly demonstrate the time scales over which dynamical phase coexistence persists.

Despite the simplicity of our model, we have shown that it is capable of describing qualitatively a 1000 year paleoecological dataset when it is given a quantity proportional to rainfall as a driving input. The key result here is that the model is able to reproduce the correlation between weather patterns and tree-grass ratios. The rough quantitative agreement between the model and data is not surprising given that these data, which are based on proxy measures instead of direct observation, are relatively crude and that only two parameters of our dynamic, nonlinear model, were tuned to achieve the fit. We note that though empirical data are frequently used in savanna modeling studies, such studies typically only consider single snapshots in time, or very short time series [START_REF] Jeltsch | Detecting process from snapshot pattern: lessons from tree spacing in the southern Kalahari[END_REF], Higgins et al. 2000, van Wijk and Rodriguez-Iturbe 2002].

While we have focused here on savannas, the concept of dynamical phase coexistence is much more general than that. Variants of the Contact Process are used in many different fields including: Physics, Ecology, Epidemiology, Sociology, etc. Furthermore, it seems likely that many of the systems to which these models are applied may feature some degree of forcing by externally varying conditions. Thus, we expect that the concept of dynamical phase coexistence will find broad applicability in a range of scientific disciplines.

  FIG.1: Schematic representation of the system dynamics. The underlying pure model has two homogeneous phases: active (woodland) and absorbing (grassland), separated by a critical point (the transition could also be discontinuous as in Fig.4). The control parameter fluctuates in time and shifts from one phase to another (droughts and rainy periods). This, combined with long-living trees, prevents the system from reaching a homogeneous steady state.

  consequently, in our model we have considered nearest-neighbor dispersal. Considering a slightly larger dispersion range (e.g., a first and second nearest-neighbor dispersal) gives rise to essentially the same type of clusters, even A c c e p t e d m a n u s c r i p t 4 if somewhat less compact. On the other hand, in the global dispersal limit clustering disappears. A larger discussion on this topic can be found, for a different but related model, in [see Calabrese et al. In press].

FIG. 2 :

 2 FIG.2: Time evolution of the tree-density in different cases. All curves except the dotted one (am = 30, σ = 6) are for am = 100, σ = 20. The dashed-dotted curve (q = 0, i.e. model without variability) exhibits damped oscillations. The solid line (q = 1) shows small variability, while the dashed and dotted curves are for q = 0.03 (intermediate variability). Notice that the one with smaller maximum age (am = 30; dotted curve) does not survive to an adverse period, while its analogous for large age (am = 100; dashed curve) does.

FIG. 3 :

 3 FIG.3: Phase diagram for the underlying pure model (q = 0) in the case of strong facilitation; observe the discontinuous transition and the hysteresis loop.

  FIG.4: Stationary tree-density ρ s vs maximum birth probability bmax, averaged over surviving realizations for the model with competition, q = 0.1, am = 100, aest = 1, and system sizes N = 40 2 (diamonds), N = 80 2 (squares) and N = 320 2 (circles).

FIG. 5 :

 5 FIG. 5: Snapshots of the system for different versions of the model with q = 0.01 and bmax = 1. (a) and (b) correspond to the same realization with aest = 4, competition, and two different times. The system fluctuates locally and globally between small and large densities, as in real savannas. (c) Facilitation, aest = 4. (d) Competition, aest = 1.

  FIG. 6: (a) Savanna mean life-time τ vs the logarithm of its size N on a log-log scale, for different values of bmax below (diamonds, bmax = 0.100 -0.104) and above (circles, bmax = 0.105-0.115) the transition point b c max 0.1045. At b c max τ grows as (ln N ) 3.68 (straight dashed line). (b) τ vs N on a log-log scale. τ diverges as N α for bmax > b cmax (circles), with the exponent α proportional to bmax (inset), indicating that the active phase is stable when N → ∞.

FIG. 7 :

 7 FIG.7: Ratio of tree to grass pollen in the Crescent Island Crater core NC93 in Kenya (dashed-dotted line) and from numerical simulations of the model (solid line), on a 100 × 100 square lattice, with age = 100. The time series of the birth probability in the simulations (dashed line) was taken to be proportional to the depth of the lake Naivasha, located next to the Island Crater, in Kenya.


Mariusz Pietruszka
email: mariusz.pietruszka@us.edu.pl

Exact analytic solutions for a global equation of plant cell growth

Keywords: growth functional, plant cells, wall extensibility

A generalization of the Lockhart equation for plant cell expansion in isotropic case is presented. The goal is to account for the temporal variation in the wall mechanical properties -in this case by making the wall extensibility a time dependent parameter. We introduce a timedifferential equation describing the plant growth process with some key biophysical aspects considered. The aim of this work was to improve prior modeling efforts by taking into account the dynamic character of the plant cell wall with characteristics reminiscent of damped (aperiodic) motion. The equations selected to encapsulate the time evolution of the wall extensibility offer a new insight into the control of cell wall expansion. We find that the solutions to the time dependent second order differential equation reproduce much of the known experimental data for long-and short-time scales. Additionally, in order to support the biomechanical approach, a new growth equation based on the action of expansin proteins is proposed. Remarkably, both methods independently converge to the same kind, sigmoid-shaped, growth description functional V (t) ∝ exp(-exp(-t)), properly describing the volumetric growth and, consequently, growth rate as its time derivative.

A c c e p t e d m a n u s c r i p t

Introduction

Plant cells are confined by a stiff, yet flexible, polymeric wall that defines cell shape and permits high turgor pressure to develop [START_REF] Cosgrove | How do plant cell walls extend?[END_REF]. In growing cells these walls must expand to enlarge and keep the wall strong enough to withstand the mechanical stress arising from internal hydrostatic pressure. The plant cell wall consists of a network of cellulose microfibrils glued together by a polysaccharide matrix (ibid.). The primary wall is secreted by the growing cell and is maintained in an extensible, plastic state during the period of cell growth, before the cell matures and the wall loses its ability to expand [START_REF] Veytsman | A model of cell wall expansion based of thermodynamics of polymer networks[END_REF]. This process, from the physical point of view, can be described by a dynamic extensibility coefficient. In this paper, three classes of cell (tissue) growth solutions are considered under time dependent extensibility coefficient. [START_REF] Lockhart | An analysis of irreversible plant cell elongation[END_REF] model deserves a prominent place among the key milestones in the field of plant growth mechanics. The model however, has many shortcomings that limit usefulness. Among these shortcomings is the formulation of uniaxial cell growth when in fact cell growth involves strain rates along three principal directions. Moreover, the mechanical anisotropy of the cell wall was ignored. This latter aspect has been considered by many authors, and quite recently a proposal was also given (see [START_REF] Pietruszka | General proof of the validity of a new tensor equation of plant growth[END_REF] and papers cited therein) concerning a new local (coordinate dependent) tensor growth equation. Yet, there is still a need for a simple global (coordinate independent), easy applicable equation of volumetric cell/organ growth keeping positive aspects of the original Lockhart equation. This paper is an effort to complement previous efforts in the area.

Lockhart equation has been known since 1965, and it is clearly not new.

However, revisiting the mathematical application of the equation can provide new insights. The Lockhart equations describe plant cell elongation, as a combination of elastic, turgor-driven extension and plastic deformation (yielding) of the cell wall. The Lockhart equation assumes that extension is driven by a constant turgor pressure due to osmotic uptake of water. The core portion of this equation has been an extensibility coefficient Φ 0 , ac-

A c c e p t e d m a n u s c r i p t

counting of visco-elastic features of cell wall. The extensibility coefficient, Φ 0 , was originally defined as a constant. This coefficient should be, however, dynamic to address the changes taking place in the walls of the growing cell of a plant. Also, numerous experiments provided evidence that the extensibility coefficient is susceptible to many environmental stimuli. These include phytohormones or growth inhibitors and temperature, and especially its pronounced change in time. From purely mathematical point of view, as well as indicated by a huge body of observational data, the extensibility coefficient magnitude should gradually decrease with time, ceasing altogether when the cell matures. In fact, the conventional treatment of the Lockhart equation (with Φ 0 = const.) delivers exclusively a divergent exponential solution for the expanding volume V , which is clearly flawed. The biological meaning of this expanding volume is that growing plant cells will extend infinitely with time in the model. Closer scrutiny of the original form of the Lockhart equation also reveals that the volumetric growth of a cell is driven exclusively by the progression of time, since the product of Φ 0 and P -Y is also a constant (see Eq. ( 1)). In other words, the cell volume expands because the time evolves. This this is indicative of a key failure of the original form of the Lockhart equation, as no inner mechanisms are responsible for growth:

ongoing biological chemistry would appear to be completely ignored.

The benchmark of success must be in better descriptions of what is observed. The experiments considering wall creep (irreversible extension) [START_REF] Cosgrove | Characterization of long-term extension of isolated cell walls from growing cucumber hypocotyls[END_REF] or where the plastic extensibility [START_REF] Cleland | The role of hormones in wall loosening and plant growth[END_REF]) experienced a kind of exponential decay with time. Unfortunately, there has been a lack of analytic description of this behavior, even at a phenomenological level, which we attempt to address herein (vide infra). Such modelling efforts are essential for description of dynamic features of the cell wall extension. Given the above arguments, a growth-equation, free of some serious drawbacks mentioned above, is desirable. In this paper a new consistent solution of these problems is presented in a form of two complementary models, namely biomechanical and biochemical which will be discussed later.

A c c e p t e d m a n u s c r i p t

Theory The Lockhart equation [START_REF] Lockhart | An analysis of irreversible plant cell elongation[END_REF] in its raw form is given by 1

V dV dt = Φ 0 (P -Y ) ( 1 )
where V denotes the volume of a growing cell. In Eq. ( 1) P and Y describe the turgor pressure and turgor threshold, respectively, constant Φ 0 stands for the extensibility coefficient. Formal solution of Eq. ( 1) gives

V = V 0 exp[Φ 0 (P -Y )t]
, where V 0 = V (t = t 0 ) is the initial volume. In a more general case, Φ 0 may be assumed as time-dependent, Φ 0 → Φ(t), and then the integration of Eq. ( 1) gives

V (t) = exp t 0 Φ(t )(P -Y )dt (2)
In the biomechanical approach, heuristically as an Ansatz having some biomechanistic underpinning, we assume that extensibility coefficient is now a time dependent function, Φ 0 → Φ(t). Time-dependent wall extensibility seems a reasonable assumption, because plant cells regulate cell wall properties by secreting enzymes, including expansins. We also assume the cell wall behaves as a damped oscillator, with wall elasticity acting as a restoring force, an inertial wall mass, and a viscotic friction due to plastic deformation of the wall. Hence we postulate that the extensibility coefficient Φ(t) satisfies the second order differential equation, able to describe elastic and viscotic (plastic) features of the cell wall, in the form

m d 2 Φ dt 2 + b dΦ dt + kΦ = f (3)
This type of equation is recognized as the time differential equation for the damped harmonic oscillator. Expanding upon the analogy with the damped harmonic oscillator, parameter m, b and k may be called inertia, viscotic friction and restoring (elastic) coefficients, respectively, and can be bound directly with the elastic/inelastic cell wall properties. Coefficient f reflecting the role of growth hormones we put momentarily as equal to zero.

Turgor pressure in Eq. ( 1) is treated as external force acting on the system.

A c c e p t e d m a n u s c r i p t

Equation Eq. ( 3) can be solved by standard procedures (strictly speaking, in non-isotropic case Φ should be a tensor). With f = 0 and Δ = b 2 -4mk > 0, the general solution has the form of aperiodic damping

Φ(t) = Φ 0 [c 1 exp(λ 1 t) + c 2 exp(λ 2 t)] ( 4
)
where

λ 1,2 = -(b ± √ Δ)/2m.
By inserting Eq. ( 4) into the modified form of the Lockhart equation

(see below) 1 V dV dt = Φ(t)(P -Y ) ( 5 )
where the cell wall extensibility coefficient depends explicitly on time, and P -Y is represented, as usual, by a constant average value, the general solution of Eq. ( 5) may be presented in the form (without the loss of generality we may put Φ 0 = 1) as:

V (t) = V 0 exp 2me -bt 2m c 1 -b + √ Δ e √ Δt 2m + c 2 -b - √ Δ e - √ Δt 2m (P -Y ) (6)
Remarkably, it turns out that Eq. ( 5) depending on the sign of the discriminant Δ one can identify three dynamic regimes (see the Appendix A for details) corresponding to various situations encountered in cell growth physiology. Some representative examples are considered in the following section.

We note that, by proposing the model Ansatz in the form of Eq. (3), we are far away from treating the cell wall literally as a damped harmonic oscillator. Rather, we use such analogy, since Eq. (3) integration produces three classes of solutions (vide infra) that we encounter in plant growth biophysics. It also turned out, that these solutions complemented by one-and two-parameter (for α-and β-expansins) biochemical models form a natural basis for simultaneous mechanical and biochemical, semi-phenomenological description of cell/organ growth physiology. This is because both approaches result in the same kind of growth functional, where the volumetric expansion V (t) is proportional to the double exponent exp(-exp(-t)) function.

A c c e p t e d m a n u s c r i p t

The biochemical model (mechanism) of cell wall loosening, due to expansin proteins action, underlying elongation growth of plant organs represents an unsolved problem and is still a matter of current debate (see e.g. [START_REF] Cosgrove | Loosening of plant cell walls by expansins[END_REF] for review). In the simplest approximation the change of expansin proteins content in the investigated system will be proportional to its concentration x in the solution and to the kinetic coefficient k 0 responsible for the interaction of endogenous expansins with the investigated sample.

Therefore, trying to tackle the problem analytically, we begin with a time differential equation for the kinetic chemical reaction in the form

dx dt = -k 0 x ( 7
)
with 0 ≤ x ≤ 1 -expansin concentration in the growing hypocotyl/coleoptile and k 0 -kinetic coefficient responsible for expansin interaction with the growing segment, leading to stress relaxation. Equation Eq. ( 7) can be easily integrated to yield

x(t) = x 0 exp(-k 0 t) ( 8 )
where x 0 = x(t = 0) describes the initial expansin concentration of the investigated system (cell or non-meristematic tissue). This equation describes the exponential degradation of an initial pulse of the wall loosening enzyme expansin. Such approach also suggests a mechanism for finite cell extension:

the cell becomes rigid as the initial pulse of expansin expires. The above argument implies equation Eq. ( 8) may be treated as a modifying factor for Φ 0 . Thus, we may substitute Φ(t) = Φ 0 x(t), which explicitly reads

Φ(t) Φ 0 = x 0 exp(-k 0 t) ( 9
)
where Φ 0 = const. is the Lockhart term (we put Φ 0 = 1 onward). Performing insertion of Eq. ( 9) into Eq. ( 5), after integration of the resulting differential

V (t) = V 0 exp - e -k 0 t x 0 (P -Y ) k 0 (10)
Above equation can be considered, similar to Eq. ( 6), as volumetric extension of a cell/tissue with V 0 denoting the initial volume.

Besides, the calculated from Eq. ( 10) half-time for stress relaxation, due to expansins action, is t 1/2 = ln 2/k 0 and can be compared with Cosgrove (1985) result (Eq. (A9), ibid.) to give k 0 = εφ. This relation, represented in our notation by k 0 = εΦ 0 , valid at least for t 1/2 , combines biochemical (k 0 ) and biomechanical (ε, Φ 0 ) aspects of plant growth, since ε stands for the Young modulus which, by definition, is a measure of the stiffness of an isotropic elastic material.

Results

Simulating plant cell growth is an important subject of research. It would be Considering the over-damped case, as in Fig. 1, we clearly see the impact of viscotic friction (b -parameter) on the shape of the function Φ (Eq. ( 14))

and consequently how it affects the growth functional V (Eq. ( 15)). This case seems to be fundamental to investigations concerning the wall creep [START_REF] Cosgrove | Characterization of long-term extension of isolated cell walls from growing cucumber hypocotyls[END_REF], see Fig. 5 therein), plastic extensibility [START_REF] Cleland | The role of hormones in wall loosening and plant growth[END_REF], Fig. 3) and growth kinetics (Edelman, 1995, Figs 1 and3) observed in

A c c e p t e d m a n u s c r i p t

various experiments. This viscotic friction parameter also appears to apply to growth retardation due to inhibitors. The latter issue is clearly visualized in Fig. 1 (right column) where with increasing coefficient b the inflection point is shifted to the higher t -values and the value of V is significantly suppressed. We believe this to be a proposal worthy of further attention, as described in the Appendix B, where the functional form of Φ(t), see Fig. 9, has been substantiated by the assumptions about the sigmoid law of great growth, stay in accord with the results presented in Figs 12.

Now, let us consider a critical case as presented in Fig. 2. This kind of behavior seems to be expected on the brink between periodic growth oscillations, as in the case of pollen tubes (vide infra) and the described above over-damped case. The resulting formula, Eq. ( 17), involved in the critical case is comparatively simple -note, it does not contain the restoring coefficient k. It looks as though, an elastic component is not present in this case and the cell wall is in a plastic state, exclusively.

Pulsatile growth appears to be the dominant mode of extension of rapidly extending tip-growing organisms. It looks as though, at this preliminary stage, the under-damped solutions to the second order differential equation, Eq. ( 3), may also find a counterpart in pulsatile growth experiments (e.g. [START_REF] Messerli | Pulsatile influxes of H + , K + and Ca 2+ lag growth pulses of Lilium longiflorum pollen tubes[END_REF] through Eq. [START_REF] Veytsman | A model of cell wall expansion based of thermodynamics of polymer networks[END_REF]. Pollen tubes from a number of plants that have been studied with adequate spatial and temporal resolution have been found to grow in an oscillatory manner. Measurements of longer pollen tubes of Lilium longiflorum that were growing in a pulsatile manner revealed a pulsatile influx of both H + and K + at the growing tip. Calcium fluxes were also measured when pulsatile growth began, the basal Ca 2+ influx decreased and a pulsatile component appeared, superimposed on the reduced basal Ca 2+ flux. The Ca 2+ oscillations have the same frequency as the growth oscillations and are nearly in phase with growth rate. It gives hope that these phenomena, by fine-tuning model parameters to receive the characteristic system oscillation spectrum (Fig. 3, mid-box), can be described in terms of time-derivative of Eq. ( 19), representing the modulated growth rate V = dV dt corresponding to Ca 2+ (K + , H + ) oscillating fluxes -see also Fig. 3 (upper plot) and compare, for instance, with by [START_REF] Messerli | Pulsatile influxes of H + , K + and Ca 2+ lag growth pulses of Lilium longiflorum pollen tubes[END_REF]. Furthermore, as we noticed, Eq. ( 5) assumes, for simplicity, a constant turgor pressure. Probably the turgor pressure will rise significantly when the wall contracts. This effect has not been included in the presented model, although it will presumably quickly damp the oscillations predicted by Eq. ( 19) (Ω-term). In order to overcome this difficulty and consider a cell like a closed volume filled with an incompressible fluid we propose a self-consistent solution (see the derivation in Appendix A). Hence, even more adequate results can be obtained with the use of an auxiliary transcendental Eq. ( 21), where more realistic system can be considered by taking into account a dampening feedback between turgor pressure and elastic wall contraction. This approach, however, demands a selfconsistent (iterative) solution as it was pointed out in Appendix A. 5) in the following form may be introduced (with the use of Eq. ( 4))

V (t) = V 0 exp t 0 c 1 e λ 1 t + c 2 e λ 2 t P (t ) -Y (t ) dt ( 11
)
where either P or Y may vary with time. This equation, among other possibilities, may also report on turgor decrease due to water deficiency, applicable in case of a previously growing plant cell that is deprived of its water resources. Thus, by accepting the adequate, specific form of P -Y term it can be solved analytically. The solution of Eq. ( 11), for the most suitable form of the exponential decay

P (t) -Y (t) = (P 0 -Y 0 ) exp(-ξt)

A c c e p t e d m a n u s c r i p t

accepted for the turgor pressure term, explicitly reads

V (t) = V 0 exp{(P 0 -Y 0 ) c 1 λ 1 -ξ e (λ 1 -ξ)t -1 + c 2 λ 2 -ξ e (λ 2 -ξ)t -1 } ( 12
)
where P 0 = P (t = 0), Y 0 = Y (t = 0) and ξ is a positive number called decay constant 1 . The reciprocal constant τ = 1/ξ is called "mean life-time"

or "scaling time" which is the time needed for the assembly (here: turgor pressure) to be reduced by a factor of e = 2.718 (Euler number). The plot of this exemplary solution is presented in Fig. 4. The remaining solutions to Eq. ( 11) are conveyed to the Appenix C.

Last but not least, Lockhart equation has been derived for elongation growth and was intended to describe a single cell elongation by assuming steady growth conditions (constant turgor). Since no constraint regarding turgor pressure is imposed onto Eq. ( 11), not only it can describe the single cell growth in varying external conditions, but can also be useful in single cell simulations. A good example is presented in the Fishman and Genard model (1998), where the growing fruit is assimilated to one big cell separated from the exterior by a composite membrane (see also [START_REF] Liu | Model-assisted analysis of tomato fruit growth in relation to carbon and water fluxes[END_REF]. Also the influence of ξ parameter, as in Eq. ( 12), onto the volumetric growth, visualized in Fig. 4, can be useful in this case.

The importance of plant growth regulators in plant tissue culture is well documented. Plant hormones (also known as phytohormones) are chemical compounds that regulate plant growth. The action of such external stimuli should also be recognized in solutions of Eq. ( 5). This can be accomplished by considering a positive f -term on the right side of Eq. ( 3). The solutions of such modified equation are shown in the explicit form in the Appendix A, as Eqs ( 22) and ( 23). An illustration of the influence of plant hormones onto the volumetric growth is presented in Fig. 5 (note, the volume V is still normalized, and does not exceed one in the case of f = 0).

It seems tempting to link both models presented in this paper. This can others, that we may obtain quantitative (numerical) values for the relevant biochemical parameter (k 0 ). More than this, since we have here considered data for the elongation growth of maize at two different temperatures, we can bind these numbers (k 0 ) directly with temperature change. The problem of temperature, in this context, is of great interest and may be subjected to rate. Temperature is a measure of the kinetic energy of a system, so higher temperature implies higher average kinetic energy of molecules and more collisions per unit time. A general rule of thumb for most (not all) chemical reactions is that the rate at which the reaction proceeds (Q10) will approximately double for each 10 0 C increase in temperature.

Once the temperature reaches a certain point, some of the chemical species may be altered (e.g., denaturing of proteins) and the chemical reaction will slow or stop. 2 Usually, an increase in temperature is accompanied by an increase in the reaction

A c c e p t e d m a n u s c r i p t

future study. Beneath we consider a preliminary result.

The environmental temperature, decisive for growth processes, can be taken into account by a well justified assumption that k 0 = k 0 (T ), where

T stands for the absolute temperature 2 . The effect of temperature on rates of reaction for a given plant species can be easily established empirically with the help of Eq. ( 10), see Fig. 8. On the other hand this outcome strongly supports the biochemical pathway resulting in the solution given by Eq. ( 10).

Discussion

Lockhart equation was conceived to explain growth for water immersed cells with negligible suction for which only turgor is required to calculate the stretching forces involved in plant growth. Our improved approach extends the idea of growth by relaxation of viscoelastic cell walls to plants subjected to suction by introducing the dynamic extensibility coefficient of the cell wall. In this sense the assumption made by introduction of Eq. ( 3) may be justified and accepted. Even though in this paper we focus primarily on the time evolution of wall samples the problem of the yield stress is not neglected. Green and coworkers (1971) provided experimental evidence on Nitella that the most labile material property in the context of Lockhart equation is the yield stress of the wall, not the extensibility. Nonetheless, it is easy to notice that all presented solutions, either of biomechanical or biochemical origin, have preserved the structure of the original Lockhart equation (also by retaining the P -Y term). Therefore, not only they can report on the wall properties, but also over the changes regarding the yielding threshold (cf. e.g. Eq. ( 11) and its exemplary solutions).

Cell expansion plays a crucial role in shaping the form and size of plants.

Expansins are plant cell proteins first discovered in studies of plant cell enlargement [START_REF] Cosgrove | Update on cell growth: Wall loosening by expansins[END_REF]. They have unique loosening effect on plant cell walls and enhance stress relaxation. Consequently, they induce a long-term, irreversible extension (creep). Even though their mechanism of action still remains enigmatic there exist temptation to link their influence with the

A c c e p t e d m a n u s c r i p t

presented model. This is made possible by invention of the analytic expressions (as Eq. ( 10)) for the wall extensibility including free parameters (x 0 , k 0 ) that can also be attributed to the action of expansins. Moreover, further equipment of our model by a gravitational force term (representing load), as a natural extension, may supposedly account for the results obtained by the experimental layout of the cell-wall extension assay [START_REF] Cosgrove | Characterization of long-term extension of isolated cell walls from growing cucumber hypocotyls[END_REF].

This in turn can easily be linked to the expansin action. This issue, concerning -among others -similar models for a more realistic case of α-and β-expansins, is presently under study.

It is worth also to emphasize that the solution of Eq. ( 5) possesses an important normalization property

V (t) ∝ lim t→+∞ exp{-exp [-t]} = 1 ( 13
)
not present in hitherto solutions. We see, we are dealing with the natural (and expected) upper bound corresponding to growth cessation when the growing cell matures. Hence, the maximum growth is treated in the model as normalized to unity. Therefore, in order to obtain significant results the control value should be always subtracted, as it is usually done in experiment. It is clear that this result extends from the over-damping to critical case as well as for all solutions of Eq. ( 10) kind.

I realize, the presented semi-phenomenological approach does not bring a new biological insight, either elucidate biological mechanisms. Phenomenology usually serves as qualitative and quantitative description of the observed phenomena. Sometimes, it leads to discovery of new solutions that can be verified experimentally or carry important hints for experimentalists or theoreticians. Also here exploited phenomenology gives such directions about changing properties of the cell wall during the process of growth. Moreover, in our case it provides interchangeable information about mechanical and biochemical aspects of expanding cell walls, also in terms of calculated numbers and empirical relations (k 0 = ln b/2) or mathematically derived rules (t 1/2 = ln 2/k 0 or k 0 = Φ 0 ). Also, easy but precise fitting procedures allow indirect (through experimental data) transformation of the following kind:

A c c e p t e d m a n u s c r i p t

(x 0 , k 0 ) ←→ (b, k) between biochemical and biomechanical "coordinates". It looks as though, with the help of the presented mathematical apparatus we can "tune" biomechanical parameters by changing the biochemical ones or, the opposite, draw conclusions about specific aspects of biochemical nature.

Also, as in the case of Fig. 7, by knowing k 0 coefficient at T = 20 0 C or at T = 30 0 C (by interpolation to the existing empirical data) we may easily predict (calculate) the plot at arbitrary temperature, below the temperature of the optimum growth (T 0 ), without time-consuming measurements. Among others, in this manner the "predicative power" of semi-phenomenological formulas can be employed.

To conclude, this paper describes a model that appears applicable to what we know empirically about the physical aspects of cell growth process. It seems clear, how realistic improvements to the modeling methods can be linked to study the biomechanical and biochemical features of cell growth. This is because of the close connection of each model with experiment through its free parameters that can be retrieved by standard interpolation of the growth data. The sparse set of parameters of both models gives also hope for a wide application in comparative studies of the growth problems encountered in plant physiological experiments.

A c c e p t e d m a n u s c r i p t

To evaluate c 1 and c 2 , the initial conditions are needed, in order to seek for these coefficients. When a damped oscillator is subjected to a damping force which is linearly dependent upon the velocity, such as viscous damping, the oscillation have exponential decay terms which depend on damping coefficient. By introducing Eq. ( 14) into Eq. ( 5) the resulting differential equation may be solved to get

V (t) = V 0 exp ⎡ ⎣ exp - b + √ Δt 2m ⎡ ⎣ 2c 1 exp √ Δ m m -b + √ Δ - 2c 2 m b + √ Δ ⎤ ⎦ (P -Y ) ⎤ ⎦ (15)
which is the explicit analytical solution for the expanding volume of the cell in the over-damped case. Here V 0 = V (t = t 0 ) describes the initial volume. The term P -Y is of the order of several tenths of MPa, and can be directly measured. Parameter b, m and Δ (k) are to be determined from experiment by interpolation of growth data by Eq. ( 15) with the help of e.g.

Levenberg-Marquardt procedure. These values then can be inserted back in Eq. ( 14) thus producing complete analytic expression for the dynamic yielding coefficient Φ(t).

Critical damping: if Δ = 0, the general solution takes a form

Φ(t) = Φ 0 (c 1 + c 2 ) exp - b 2m t ( 16
)
and c 1 and c 2 are to be obtained from initial conditions. By inserting Eq. ( 16) into Eq. ( 5) the evaluation of the differential equation results in

V (t) = V 0 exp - 2m b (c 1 + c 2 ) exp - b 2m t (P -Y ) (17)
Again, this is the explicit analytical solution for the expanding volume of the cell in the critical (maximum damping) case. Here V 0 describes the initial volume. The only parameter b and m are to be determined from experiment by interpolation of growth data with Eq. ( 17). These values then can be plugged in Eq. ( 16) to receive analytic expression for the dynamic yielding coefficient.

A c c e p t e d m a n u s c r i p t

Under-damping (hereafter referred also to as cell growth (rate) oscillations):

if Δ < 0, the solution reads

Φ(t) = Φ 0 exp -b 2m t A cos |Δ| 2m t + B sin |Δ| 2m t ( 18
)
where from the given initial conditions, the coefficients A and B can also be evaluated numerically. By inserting Eq. ( 18) into Eq. ( 5) the integration of the resulting differential equation yields

V (t) = V 0 exp m exp -bt 2m b 2 -2mk (P -Y ) [Ω] (19)
where

[Ω] = -(Ab + B √ Δ) cos t √ Δ 2m + (A √ Δ -Bb) sin t √ Δ 2m
Eq. ( 19) is the explicit analytical solution for the expanding volume of the cell in the oscillatory case. It exploits the sigmoid growth double exponent in the form of exp(-exp(-t)) in a superposition with a periodic function (Ω -pulse); b, m, k are free parameters to be determined from experiment by interpolation. Constants A and B contribute to the superimposed pulse amplitude whereas √ Δ 2m = ω = 2πν -represents its frequency ν in Hz, respectively. These values can be inserted in Eq. ( 18) to get an expression for the dynamic yielding coefficient in this case.

A more subtle approach requires introduction of a dampening feedback between turgor pressure and elastic wall contraction. In order to bring this feedback between turgor pressure and elastic wall contraction into Eq. ( 5)

we utilize the state equation (we accommodate argumentation as given by [START_REF] Stanley | Introduction to phase transitions and critical phenomena[END_REF] and restrict to the first approximation in the virial theorem).

Thus, we may assume P = γT /V where T is the absolute temperature (in Kelvin scale). The latter assumption needs a comment. There exist a few V is generally interpreted as cell volume but in the context of numerical recalculation of experimental data in light of our model we may also refer

V to the entire organ, the coleoptile or hypocotyl. The pressure of inner cell solution is proportional to T and this is a natural consequence of the molecular energy increase with temperature3 . The reciprocal proportionality of P to V we introduce as in the first approximation, which gives the linear solution V = V (t) of Eq. ( 5). This assumption goes along with the fact that the experiments are usually performed in the linear range of the sigmoid growth curve. It is, however, evident that in more precise calculations one should solve Eq. ( 5) for subsequent reciprocal powers of the volume V by the iteration method. Such extension would be even more adequate, nonetheless also more complicated. However, that high accuracy in case of biological experiments where we deal with relatively high statistical error is superfluous.

Giving the above arguments together, we receive the following equation

1 V dV dt = Φ(t) γT V -Y (20)
instead of Eq. ( 5), where T is the absolute temperature and γ is a constant.

Returning to the problem of cell oscillations we may also change the right side of Eq. ( 19) in analogous fashion, replacing P by γT /V , to receive the self-consistent transcendental equation for the fluctuating volume V : (21) may be solved selfconsistently by iterations. It looks as though, this is the simplest way the dampening feedback between turgor pressure and elastic wall expansion (contraction) may be taken into account.

V (t) = V 0 exp m exp -bt 2m b 2 -2mk γT V -Y [Ω] (21)
Stimulator present One of the most important solutions, taking into account the role of growth regulators (phytohormones), can be obtained by adding a non-vanishing exerting force f on the right side of Eq. ( 3). Then, the solution of Eq. ( 3) reads

Φ(t) = Φ 0 c 1 exp - b - √ Δ 2m t + c 2 exp - b + √ Δ 2m t + f k (22)
and corresponds to the volume extension in the form

V (t) = V 0 exp exp - b + √ Δt 2m [Υ] + f t k (P -Y ) (23)
where

[Υ] = ⎡ ⎣ 2c 2 exp √ Δt m -b + √ Δ - 2c 1 m b + √ Δ ⎤ ⎦

Appendix B

Independent estimation of the dynamic yielding coefficient

Based on the Lockhart model and the assumption that the cell volume variation follows a sigmoid curve a formula for the extensibility coefficient can be constructed. Thus, one can explain, in a semi-phenomenological manner, the shape of the cell wall extensibility experimental plots encountered.

In plant physiology growth is described by the law of great growth Fogg Recalling the Lockhart equation Eq. ( 1) and assuming Φ 0 as timedependent Φ 0 → Φ(t), integration of Eq. ( 1) yields

V (t) = V 0 exp t 0 Φ(t )(P -Y )dt (24)
On the other hand, assuming sigmoid character of growth, one gets

V (t) = 1 2 V 0 1 + tanh t t 0 -1 (25)
By comparing Eqs (24) and (25), assuming the hydrostatic pressure constant P -Y =const=c, calculating log of both sides and taking the time-derivative, we arrive at an analytic expression for the dynamic yielding coefficient in the form

Φ(t) = sec 2 t t 0 -1 ct 0 1 + tanh t t 0 -1 (26)
In the above equation t 0 denotes the inflection point of the elongation/growth curve that corresponds to the time of the maximum extension rate. The dimension of Φ follows: [Φ]=(MPa) -1 (time unit) -1 , consistently with the Lockhart equation. Providing the application of a proper time-scale (minutes, hours or days) to a particular growth problem, the above formula can be interpreted as an approximated expression for the dynamic yielding coefficient Φ (see also Fig. 9).

A c c e p t e d m a n u s c r i p t

Appendix C

Special solutions to Eq. ( 11) -Exact results

Among others, we may consider the following special cases.

1. P (t) -Y (t) = P 0 exp(-ξ 1 t) -Y 0 exp(-ξ 2 t), with separate decay constants ξ 1 and ξ 2 by P and Y processes, respectively. 19) resembling tip-measured pollen tube influx of Ca 2+ -a periodic signal with random noise added. The tubes show pulsatile growth (lower plot) and pulsatile Ca 2+ infux (upper plot). The Ca 2+ oscillations have the same frequency ω as growth oscillations and are nearly in phase with growth; that is, the peaks of tip-localized Ca 2+ nearly coincide with growth rate. The calculated spectrum via Fourier transform (mid-plot) shows a very strong peak at 30, the calculated "frequency" of the calcium influx originating the periodic signal and oscillatory growth response. There is a second peak, which is essentially a consequence of aliasing. 10), to the original data (*, +) taken from [START_REF] Lewicka | Theoretical search for the growth temperature relationship in plants[END_REF], Fig. 2. The expansin pool is assumed to be at the same level for the initial time, corresponding to x0 = 1. The resulting kinetic coefficient k0 leading to wall stress relaxation and subsequent volumetric growth equals: k0 = 0.252611 (smaller) for T = 20 0 C and k0 = 0.418721 (greater) for T ≈ 30 0 C , as expected. The determination coefficient R 2 for the fitting procedure exceeds 0.99 (R 2 allows us to determine how certain one can be in making predictions from a considered model: here a differential equation V -1 V (t) = x0 exp(-k0t)(P -Y ) leading to Eq. ( 10)). Left scale in [mm] -we take the coleoptile cross-section equal to 1 mm 2 . The procedure is valid for all temperatures up to the optimal growth temperature.

V (t) = e c 1 ( -1+e t (λ 1 -ξ 1 ) ) P 0 λ 1 -ξ 1 + c 2 ( -1+e t (λ 2 -ξ 1 ) ) P 0 λ 2 -ξ 1 - c 1 ( -1+e t (λ 1 -ξ 2 ) ) Y 0 λ 1 -ξ 2 - c 2 ( -1+e t (λ 2 -ξ 2 ) ) Y 0 λ 2 -ξ 2 V 0 (27) 2. P (t) -Y (t) = P 0 -Y 0 exp(-ξ 2 t), P 0 = P (t = 0) = const. V (t) = e c 1 ( -1+e λ 1 t ) P 0 λ 1 + c 2 ( -1+e λ 2 t ) P 0 λ 2 - c 1 ( -1+e t (λ 1 -ξ 2 ) ) Y 0 λ 1 -ξ 2 - c 2 ( -1+e t (λ 2 -ξ 2 ) ) Y 0 λ 2 -ξ 2 V 0 (28) 3. P (t) -Y (t) = P 0 exp(-ξ 1 t) -Y 0 , Y 0 = Y (t = 0) = const. V (t) = e c 1 ( -1+e t (λ 1 -ξ 1 ) ) P 0 λ 1 -ξ 1 + c 2 ( -1+e t (λ 2 -ξ 1 ) ) P 0 λ 2 -ξ 1 - c 1 ( -1+e λ 1 t ) Y 0 λ 1 - c 2 ( -1+e λ 2 t ) Y 0 λ 2 V 0 (29)

A c c e p t e d m a n u s c r i p t equation, we arrive at the formula

  especially important when we need to further improve the plant cell growth properties in the future. Quite some research effort has been devoted in this research area. This article introduces an extension of the Lockhart equation that greatly expanded the predictive capacities of previous equation. Therefore, we are motivated first to find the basic semi-phenomenological solutions of growth equation that provide a better description of the time dependent extensibility coefficient. The interplay of parameters in our model results in three kinds of solutions (ready to use formulae are conveyed to the Appendix A). First exemplary results, due to this classification, are presented in Figs 1-3.

  Figs 4-5, 7 in paper A c c e p t e d m a n u s c r i p t

  Another question arises concerning constant P -Y term. The assumption that P and Y are unchanging (P -Y = const.) is valid only under steady growth conditions. What are the consequences of this assumption to apply the equation at a long time scale under fluctuating conditions (say, over few weeks)? For instance Lechaudel et al. (2007) reported a nonconstant threshold pressure Y during fruit growth. To answer this question the general solution of Eq. (

A c c e p t e d m a n u s c r i p t be accomplished with

  Fig.1a-cand then interpolated with the help of Eq.[START_REF] Green | Metabolic and physical control of cell elongation rate[END_REF]. As a result of the fitting procedure (x 0 , k 0 ) pairs are obtained and new plots (Fig.6) are created which reproduce the original data from Fig.1. This easy fitting exercise shows a good agreement between both sets of growth curves (compare plots in Fig.1a-c and Fig.6), and consequently between its biomechanical and biochemical representations. By setting x 0 as fixed (x 0 = 1) we may repeat calculations to receive k 0 = 0.602, 0.369 and 0.309 corresponding to b = 2.1, 3 and 3.5, respectively. A closer scrutiny reveals (see the caption) that the greater k 0 coefficient the lower b value, and the opposite, as it should be expected. Another extrapolation based on above presented numbers approximately gives k 0 = ln b/2, binding mechanical (b) and biochemical (k 0 ) model quantities together in one empirical formula. Needless to say, exactly the same routine can be performed with data taken directly from any real experiment, giving the opportunity to interconnect biochemical and biophysical aspects (parameters) of volumetric plant cell growth (elongation) for a given species. Similar procedures can be conducted, depending on the experimental context, with the remaining equations presented in this paper, thus strengthening the predicative power of the presented model. Also, in order to support our analysis let us take some real, experimental growth data from Lewicka and Pietruszka (2006) -see Fig.2therein. The result of calculations is given in Fig.7(see the caption). It shows, among

A c c e p t e d m a n u s c r i p t able to propose

  empirical data (e.g.[START_REF] Proseus | Turgor, temperature, and the growth of plant cells: using Chara corallina, as a model system[END_REF] revealing weak turgor threshold Y dependence on temperature T . Despite this sparse information we are not an exact model function for Y (T ). Nevertheless, in zeroth approximation we may surely assume Y (T ) = const. Based on this work (ibid.) we acknowledge that the Y (T ) dependence is much smaller than P (T ) and as such may be neglected (the assumption Y = const. still holds).

A c c e p t e d m a n u s c r i p t

  Above equation can be solved by simple iterations (or other methods: bisec-tion, regula falsi, secant, Newton-Raphson). The temperature T is, of course constant during experiment, therefore we may put P = const/V in Eq. (20), and consequently also in Eq. (21). Then, the transcendental equation Eq.

(A c c e p t e d m a n u s c r i p t

  1975). There is also common agreement that because the sigmoid function has a non-negative first derivative in positive domain and exactly one inflection point, this functional form properly reproduces the large-scale evo-lution in time. In fact, growth of any plant organ can be split into three basic phases: the initial phase of slow growth, the intense growth phase and, eventually, the final phase of slow growth. Such regularity can be represented by a sigmoid curve that characterizes the time evolution of individual cell growth and the growth of plant organs as a whole.

Figure 1 :

 1 Figure1: Dynamic cell wall extensibility Φ, as defined by Eq. (14), and the corresponding volumetric growth V , Eq. (15), as a function of time t. Over-damping case. The parameters are:P -Y = 1, m = 1, k = 1, b = 2.1 (a and a'), b = 3 (b and b'), b = 3.5 (c and c'). All figures are described in terms of suitable units for the considered time-scale; also we adopt such units that Φ0 = 1. We put c1 = c2 = 1 in Figs 1-2, 4-5 and V0 = 1 in Figs 1-4. It is easy to notice how slightly changing wall mechanical properties (bparameter) drastically change the shape of the growth functional V .

A c c e p t e d m a n u s c r i p t 2 Figure 2 :

 22 Figure 2: Dynamic cell wall extensibility Φ, as defined by Eq. (16), and the corresponding volumetric growth V , Eq. (17), as a function of time t. Critical damping case. The parameters are: P -Y = 1, m = 1, b = 2.

Figure 3 :

 3 Figure 3: Numerical simulation (A = B = 1, Δ < 0) based upon Eq. (19) resembling

Figure 4 :

 4 Figure 4: Volumetric growth V as defined by Eq. (12) as a function of time t. The parameters are: m = 1, k = 1, b = 3, ξ = 0 (upper curve), ξ = 0.05 (lower curve)applicable in case of a previously growing plant cell that is deprived of its water resources due to water deficiency (slightly exaggerated to highlight the effect).

Figure 5 :AFigure 7 :

 57 Figure 5: Volumetric growth V (upper plot), in presence of growth regulators, as defined by Eq. (23), as a function of time t. Net growth V -V0 (lower plot), where V0 corresponds to f = 0 -control, as an example of linear reaction due to the use of growth stimulator. The parameters are: P -Y = 0.8, m = 1, k = 1, b = 2.1; (a) f = 0 (no stimulator), (b) f = 0.005 and (c) f = 0.01

Equation Eq. (12) will be finite for long times if ξ ≥ λ1, λ2. Then we receive:limt→∞ V (t) = exp{(P0 -Y0) c 1 ξ-λ 1 + c

ξ-λ 2 }.

This statement stays also in agreement with the state equation in the general form ofP = wρ = w k B TV , where ρ is the energy density and w is a constant; kB is the Boltzmann constant.

Acknowledgment

I kindly thank Dr P. Gusin, Institute of Physics, University of Silesia for critical reading of the manuscript.

Appendix A

Special solutions -exact results

Based on the sign of the discriminant Δ solutions to Eq. ( 3) may be divided into three categories.

Over-damping case: if Δ > 0, the solution is written explicitly as (input data: see the caption to Fig. 7). The optimum growth temperature for maize approximately equals T0 = 300 K. The kinetic coefficient k0 changes in the interval of 20-30 0 C about two times. This is in agreement with "Q10 law" stating the rate at which the reaction proceeds will approximately double for each 10 0 C increase in temperature.

The procedure is valid in ascending part of growth/elongation vs temperature curve (below T0).