Plant behaviour: a mathematical approach for understanding intra-plant communication

Plants are far from being passive organisms being able to exhibit complex behaviours in response to environmental stimuli. How these stimuli are combined, triggered and managed is still an open and complex issue in biology. Mathematical models have helped in understanding some of the pieces in the complexity of intra-plant communication, but a larger and brighter view, setting together multiple key processes, is still missing. This paper proposes a fully coupled system of nonlinear, non-autonomous, discontinuous, ordinary differential equations to describe with accuracy the adapting behaviour and growth of a single plant, by deeply analysing the main stimuli affecting plant behaviour. The proposed model was developed, and here sustained, with the knowledge at the state of the art; and validated with a comparison among numerical results and a wide number of biological data collected from the literature, demonstrating its robustness and reliability. From the proposed analysis it is also shown an emerging self-optimisation of internal resources and feedback stimuli, without the need for defining an optimisation function for the wellness of the plant. The model is ultimately able to highlight the stimulus-signal of the intra-communication in plant, and it can be expanded and adopted as useful tool at the crossroads of disciplines as mathematics, robotics, biology, for instance for validation of biological hypothesis, translation of biological principles into control strategies or resolution of combinatorial problems. Author summary Plants are able to adapt themselves to a wide range of conditions. The complexity and efficiency of this behaviour suggest a richness of signals exchanged that allowed scientists to talk about communication. Understanding this network of stimuli is an important challenge to increase biological knowledge as well as to adapt plant mechanisms to different research fields. To address this quest we developed a mathematical model able to describe what cues are fundamental in plant communication and how they interact each other.

1 Introduction morphology in response to the environmental changes is referred by biologists as plant 5 behaviour [2], or plant plasticity [3]. 6 In a more recent review, [4] redefined the definition of plant plasticity adding to the 7 previous environment-induced stimuli, the ability of plants to live according to 8 purposeful behaviours developed on the base of a sort of memory for some fitness 9 optimisation. Moreover, according to the author, the large number of external stimuli 10 and internal targets that the plant has to manage could be the reason of not placing 11 critical functions in a single (or two) organ, like animals' heart or brain. 12 Hence, plants are able to trigger a network of different cues that control molecular 13 signalling in different tissues, organs and the whole plant [5]. Decoupling and clearly 14 understanding such signals are challenging tasks but of fundamental importance for 15 plant behaviour analysis. According to the studies that better analysed plant 16 behaviour [1] , intra-plant signals can be categorised into three main groups: I) electric 17 signals, typically used by plants to stimulate a fast adaptation, as required for instance 18 in case of some nutrient deficiency or wounds [6,7]; II) oscillatory signals, that are 19 driven by the circadian clock, and direct the growth [8,9]; and III) hormons, like auxin, 20 citokin, giberellin, which are typically messengers of stress conditions [10,11]. 21 Towards the understanding of plant behaviour, several mathematical models have 22 been proposed with the aim of analysing selected signalling pathways and internal 23 dynamics. For example, in [12], a coupled PDE-ODE (Partial and Ordinary Differential 24 Equation respectively) system describes the cross-talk between the brassinosteroid and 25 gibberellin signalling pathways; while, other authors modelled the water uptake by roots 26 and the associated morphology [13,14]. The fluid dynamics inside the xylem and 27 phloem has been also of interest [15,16], as well as, the photosynthesis [17][18][19][20] and the 28 interaction between different plants [21,22]. However, plant behaviour is determined by 29 a very complex assembling of very different signals and to date none of the proposed 30 models is able to describe such complexity and clearly discriminate 31 stimulus-signal-behaviour chain in view of an internal communication network. 32 A previous tentative of analysing the mechanisms behind intra-plant communication 33 has been done in [23], with the aim to identify and exploit basic biological principles for 34 the analysis of plant roots behaviour. The analysis was carried in simulated 35 environment and on robotic roots by translating the identified biological rules into 36 algorithmic solutions. Particularly, on each root individually, the uptake kinetic and a 37 local memory were used to adjust the local priority of nutrients, while the resources 38 collected as a whole were shared according to the source-sink principle. Even if each 39 root is acting independently from the others, just considering local information, results 40 shown the emergence of a self-organising behaviour in plant roots aim at optimising the 41 internal equilibrium among nutrients at the whole plant level, with no need for a central 42 coordination or an optimisation function. However, in that work just a small part of the 43 plant complexity was considered, for instance, it completely excluded the analysis of the 44 above ground organs and photosynthesis-related dynamics. 45 In this paper, we do not aspire having a complete description of the plant 46 complexity, yet we want to identify the main cues influencing the growth of a plant with 47 the aim of investigating on the processes playing a role in the intra-communication for 48 plant growth decisions. We propose and explain here a system of ODEs that, differently 49 from the state of the art models, take into account the entire sequence of processes from 50 nutrients uptake, photosynthesis and energy consumption and redistribution. 51 In particular, according to [24], biomass is an important estimator of plant fitness 52 but it cannot be studied alone. For this reason, our model focuses on both above and 53 below ground biomass, along with nutrients and sugar levels. In this way, it is possible 54 to study independently leaves and roots adaptations to external and internal conditions. 55 Moreover, unlike the most of existing models, it is possible to describe plant behaviour 56 when a more complete interaction among nutrients is taken into account. 57 We propose here a system of ODEs to outline the main cues influencing the growth 58 of a plant. In section 2 the model is described in details, emphasising the well proved 59 biological assumptions used to build up the equations. Section 3 presents several results 60 of numerical simulations and their comparison with the biological data obtained from 61 the literature, for validation of the model. Finally, in section 4 are discussed the key 62 elements of our model, its functionalities and the improvements with respect to the 63 state of the art, while section 5 presents conclusive remarks. 64 2 Methods 65 2.1 Description of the mathematical model 66 Through the photosynthesis, plants produce glucose from carbon dioxide and water, 67 using the energy from sunlight. To avoid starvation in the night, they should carefully 68 manage the resources acquired and produced during the day. Also, a plant is composed 69 by organs with different tasks. For instance, roots are the organs of the plant furnishing 70 water and nutrients to the above ground part, receiving from it the photosynthetic 71 products, particularly in the form of sugars, needed for growth and maintenance. While 72 the leaves are the energy producers and starch stores used during nights. According to 73 the resources availability and the energy produced, the biomass is differently distributed. 74 Nutrients are transported toward root surface mainly by water flow in soil. 75 Nevertheless, nutrient uptake is controlled by spatial distribution of roots. In fact, as 76 explained in [25], root system can develop to easily access to immobile nutrients (e.g. 77 phosphorus or ammonium), or locally adapt their exploitation power for a specific 78 nutrient. 79 According to [26], there are 17 most important nutrients for plants that can be 80 distinguished in two main categories: macronutrients (usually required in high 81 quantities) and micronutrients (which are required in smaller quantity). 82 As reviewed in [27], the nutrients content in soil should fall in specific ranges for the 83 plant not to suffer from toxicity or deficiency effects. Moreover, the same authors show 84 how an optimal ratio between nutrient contents in plant is necessary and how an 85 imbalance of such ratio can actively affect the uptake of nutrients. 86 Due to the complex network of interactions between nutrients, in both soil and plant, 87 usually only one nutrient [17,20,28], or the single water availability [14,21,29], have 88 been studied. 89 Some authors have studied an interaction between nutrients but avoiding a whole 90 plant interaction like photosynthesis or sugar-signalling [30].

91
In our model, in order to study a more complete network of signals, over 92 sugar-signalling, we developed a root process dynamic model considering the presence of 93 two nutrients in the soil, and assuming a non-limiting and uniform distribution of water 94 in soil and of CO 2 in air.

95
For our study, we selected nitrogen (n) and phosphorus (p), being both fundamental 96 nutrients for growth in plants. 97 Nitrogen is the most abundant nutrient in plants, essential constituent of protein 98 and chlorophyll. Moreover it promotes root growth, improves fruit quality, enhances the 99 growth of leafy vegetables, increases protein content of fodder, encourages the uptake 100 and utilisation of other nutrients [31].
The proposed model takes in consideration all these fundamental working principles, 106 by describing the dynamics of each actor: photosynthesis, starch, sucrose, nutrients 107 uptake and management, costs of maintenance and growth, and biomass allocation.

108
For the sake of simplicity, the time dependence of functions and variables is implied 109 and it will not be reported in the following model equations. 110

111
We described the photosynthesis ([ µmolC6 gF W ·h ]) process with a time-dependent function describing the content of sugar produced by each gram of photosynthetically active leaf biomass: where:  [34] suggests a more specific photosynthesis limiting process, which is indirectly linked to the need of the plant of avoiding the production of starch and sugar in excess. Accordingly to [35] and [20], we propose a limiting function depending on the maximum starch that the plant can consume in nightly hours. The C function can be written as: where: The parameter λ c describes the strength of photosynthesis limitation, a is the 117 starch, a max the maximum starch that can be consumed at the maximum rate 118 τ max as during the night and f the photo-period. The sugar limitation will affect the 119 starch partition (γ) dynamics, that will be explained later (Section 2.3.2).

120
• p max h [ µmolC6 gF W ·h ] is the maximum rate of photosynthesis. As for in [36], it is fixed to 121 12.7 µmolC6 gF W ·h .

122
• n t and p t represent nitrogen and phosphorus saturating threshold, respectively.

123
Both these nutrients play a role in the photosynthesis process (see [37] for an 124 explanation of the role of nitrogen, while [38] reviews phosphorus influence on 125 photosynthesis). Hence, photosynthesis is strictly correlated with nutrient 126 availability in a saturating way. In fact, according to the law of minimum [39], the 127 most limiting nutrient content will mainly affect the photosynthesis. By [17], the 128 rate of carbon production obtained from nitrogen content is estimated as: Hence, the minimum nitrogen content required to sustain the maximum rate of photosynthesis is: The saturating function n t can be computed as: Not having found similar information for phosphorus, we estimated the saturating function of this nutrient by assuming the consumption of phosphorus as one-tenth of the nitrogen's one, according to the optimal ratio among these two nutrients, which has been observed to be about 10 (as reviewed in [23]). Therefore, we approximated: In particular, if n < n ph or p < p ph , then it should be n t = 0 or p t = 0 respectively just 132 because there are not enough resources to start the chemical processes.

133
It should be noted that Eq (1) is a simplification of the more general formula: which takes into account the aboveground biomass (b l ) and the maximum

139
During the day, some photosynthate is stored as starch to be degraded during night to sustain the nocturnal metabolism. Independently by the length of the night, the starch is degraded into sucrose in a nearly linear manner (at a rate τ as ), such that almost all of the starch is used by dawn [40]. We described starch and sucrose dynamics with the equations: A very similar dynamic can be found in [41] and [42]. Nevertheless, some important sugar-signalling and leaf competition (as described in previous section); 144 • The starch partition coefficient γ is not a constant, neither it is estimated 145 assuming sucrose homoeostasis (as in the model of [43]). It provides instead plant 146 adaptation as a function of time (its dynamic is described in Section 2.3.2);

147
• The starch degradation rate τ as is simpler respect to the definition provided 148 in [41] and [42], and it does not depend on a subjective dusk [44]. Even 149 considering the simplification adopted in this paper for starch degradation, its 150 dynamic remains correct (see file S1 in supporting information for demonstration); 151 • Just like in [20], both uptake (r u m ) and transport (r t m ) processes are sucrose 152 consuming, over the usual maintenance respiration (r m );

153
• The rate of sucrose sent into the phloem intended for growth (ηr g ) is affected by 154 nocturnal efficiency [45]. It is well established that maintenance and growth respiration occur both during day and night. The sucrose to sustain nightly metabolism is provided by the degradation ([ µmolC6 gF W ·h ]) of starch stored in the morning [40], and here defined as: where:

157
• τ max as = 6 µmolC6 gF W ·h is the maximum degradation rate estimated by the results

159
• s max = 2 µmolC6 gF W is the maximum sucrose content in leaves. It is estimated by 160 experimental values in [46]. A negative feedback on starch degradation due to 161 high sucrose levels is proposed also in [47]. Hence, the introduction of a maximum 162 threshold is justified by these evidences. On the other hand, a minimum threshold 163 should be expected to trigger sugar production or starch degradation in case of 164 starvation. Like in [36], s min is estimated as 1.3 µmolC6 gF W ;

165
• a | dusk [ µmolC6 gF W ] is the starch content at the beginning of the night. The use of this 166 parameter makes us highlight some differences with respect to the model 167 presented in [41], where: I) the function proposed to describe starch degradation 168 is independent from this parameter and has a peak at dawn; II) starch 169 degradation is assumed to exist also during the light period; and III) the non 170 linear discontinuous function proposed takes into account the surface of starch 171 granule. In the supporting material, we integrated their function in our model and 172 we compared the two approaches;

173
• a min = 0.15 µmolC6 gF W is the minimum amount of starch that plant has to ensure. It 174 is estimated by [41]. 175 Hence, during night (L = 0), the starch is degraded with a rate that depends both on 176 the time t with respect to the next expected dawn (circadian clock), the sucrose 177 signalling and the constant rate necessary to consume almost all of the starch by dawn. 178 This strategy holds unless the required rate is greater than the maximum one.

Starch and sucrose partitioning
180 [47] explains how rising levels of sucrose stimulate starch synthesis, while daily sucrose starvation decreases starch accumulation. Nightly sucrose starvation, on the other hand, promotes starch production [20]. From these observations, we defined the dynamic of starch-sucrose partitioning as:

181
Sucrose transport is essential to allow growth and maintenance of non 182 photosynthetically active tissues like roots, shaded and/or young leaves, and other 183 storage organs [48]. Actually, the sucrose is loaded into phloem generating an osmotic 184 pressure which drives both water and sucrose from source tissues to sinks [49].

185
In particular, testing Münch hypothesis, [50] shows a sink priority concept affecting how sucrose loaded is allocated between sinks. In our model, sucrose is divided between above and below ground biomass. The partitioning coefficient f r of sucrose allocated to roots is described by the following equation: where f n = 10 10+ n p is a stoichiometry signal with respect to the optimal nutrient ratio.

186
The frequency parameter fixed to 1 1 h is omitted in the previous equation.

187
According to [51], the less nutrients are, the more root priority is and vice versa. In 188 particular, [52] stresses that sink priority is affected by the most limiting resource, 189 while [53] and [54] highlight how plants are able to adapt their allocation strategies 190 according to internal and external stimuli. Hence, while functions a n and a p describe 191 nutrients demand, the value of s with respect to the threshold s min describes a 192 limitation due to leaves.

193
As result, the function f r balances the nutrients demand with respect to the most 194 limiting resource, and describes the sucrose allocation between leaves and roots.

195
Moreover, being a central brain missing in plant, there is not central control on 196 carbohydrates allocation, thus we hypothesis a different nutrient feedback driving 197 resources allocation within roots (the same hypothesis should be valid for leaves, but we 198 are not interested in modelling spatial distribution of leaves at this stage of the work). 199 According to [55], [51] and [56], the innate soil heterogeneity affects root system morphology promoting root proliferation in nutrient-rich soils. In our model the choice is between two different soil zones. Given sucrose allocated to roots, the proportion of resources sustaining growth in the top soil is given by the following weighted average Being soil zones vertically distributed, the binary function χ 1 is necessary to express if 200 the root biomass has reached the second (deeper) soil zone. Actually, r is a parameter describing the critical biomass to reach before having access to 202 the deeper zone.

204
Considering the above observations on biomass allocation, we can describe the growth dynamics of leaves and roots with the equations: In Eq (16), (17) and (18), λ sb is a conversion parameter from sucrose to biomass. The 205 growth, according to the specific tissue priority (f r and e 1 ), will be proportional to the 206 whole sucrose exported by photosynthetically active leaf biomass. To study the root system adaptation in a changing environment, we assumed the soil 209 divided into two zones: the topsoil and the subsoil. This choice is justified by the 210 different mobility of phosphorus and nitrogen (the selected nutrients for our study).

211
The former is immobile and usually a greater content can be found in topsoil [57].

212
Nitrogen, on the other hand, is quite mobile and easily transported by percolating water 213 in subsoil [58,59]. 214 We described the plant dynamics of nitrogen n [ µmolN gF W ] and phosphorus p [ µmolP gF W ] content as: In [17], the cost of assimilating n into biomass is estimated as As already previously justified (Section 2.2), the p cost of assimilating C 6 can be approximated using the optimal ratio value ( [23]) as c sp = csn 10 . In particular, metabolism, transport and growth will be reduced if the nutrient content is not sufficient to sustain the costs (c sn and c sp ). Then: The second terms of Eq (19) and (20) are the n and p costs owed to the use of sucrose 220 for maintenance and growth. Also, the starch requires a maintenance respiration (and 221 then a cost in nutrients) due to the production of storage structures [60].

222
Note that s and a represent the quantity of moles of C 6 for each gram of 223 photosynthetically active leaf biomass. Then, the term is necessary to 224 consider the actual plant sugar content, being B the full plant biomass. 225 Finally, the last terms of Eq (19) and (20) express the nutrient costs due to 226 photosynthesis. Knowing the nutrient cost for the maximum rate, n ph (Eq (5)) and p ph 227 (Eq (7)), the actual consume will be proportional to the intensity of photosynthesis.

228
The frequency parameter λ f [ 1 h ] will be estimated later (Section 2.9). Having nutrients located in two discrete soil zones, we described the potential uptake with two terms (labelled in subscript with index 1 -topsoil -and 2 -subsoil), each multiplied by a factor that takes into account the whole biomass. For nitrogen (and similarly for phosphorus) we defined: where b r1 is the root biomass in topsoil, b r2 the root biomass in subsoil. b max r1 and b max r2 231 are the maximum root biomass allowed in each soil zone. As for b max l (in Eq (9)), we fix 232 to infinite their value, assuming no root competition. a n is an internal feedback 233 explained later.

234
A limiting factor [−] should be taken into account. Its role is to decrease the forage if the hourly potential uptake u pot n needs excessive sugar consumption with respect to the available one. Thus: u n1 and u n2 [ µmolN gF W ·h ] are the nutrient rates of uptake, while n c will be estimated later. It is well established that the uptake follows a Michaelis-Menten kinetics [61,62]: constant and n s1 [ µmolN cm 3 ] is a time depending function describing the nutrient content in 236 the soil zone. Similar equations are defined for u n2 and phosphorus in top and subsoil. 237 Several papers have shown how Michaelis Menten parameters are not fixed [63][64][65]. 238 Plants, in fact, are able to modify these parameters to be more or less affine with the 239 specific nutrient. This adaptive behaviour is experimentally measured changing nutrient 240 soil content [66]. 241 However, [55] notes how changes in kinetics parameters could be due to plant 242 internal status instead of environmental reasons. In particular, [55] reports how root 243 proliferation and uptake require varying energies, and how the balance between costs 244 and gains could explains the plant strategy.

245
In our model we tried to describe this behaviour. For this reason, from [67] we estimated kinetic parameters for nitrate, and from [63] we approximates phosphorus parameters: They represent the maximum values for the couple (I max , k). To simulate the nutrient affinity adaptation we introduced in Eq (25) the internal control a n [−]. An analogous control a p is inserted for phosphorus dynamics Both a p and a n are positive time dependent functions bounded by 1. Their dynamics are described below (again the frequency parameter fixed to 1 1 h is omitted in both the following equations): A similar formulation holds for the phosphorus control function a p : C n and C p represent the nitrogen and phosphorus costs due to the actual strategy. They 246 are the absolute value of the sum of the second and third terms in Eq (19) and (20). 247 n max and p max take into account the memory of a plant. In fact, according to [18], plants manage its foraging strategies in order to avoid both scarcity and excess of nutrients in stores. In particular, the optimal nutrient status is defined by how many days plant can survive maintaining the same rate of growth if a specific nutrient is not available anymore. [18] estimated for herbaceous plants (like Arabidopsis) a memory of 4 days. We fixed D = 4 [−]. Hence, the minimum amount of nitrogen content ([ µmolN gF W ]) required to sustain photosynthesis and respiration for one day at the maximum sugar content is: Then, we estimated n max = Dn min and p max = n max 10 . The term: (1 − u n u n + C n ) n max n max + n describes the active incentive for the uptake rate according to the internal status and the cost-gain balance as reported in [55]. On the other hand, the last terms in Eq (29) and (30) represent the active limitation to the uptake due to the same concepts. Moreover, as emphasised by [68], over an active control of nutrient affinity, the uptake rate is passively increased by photosynthesis and leaves transpiration. The term: (32) describes this behaviour. 248

249
All the sucrose s produced by photosynthesis and starch degradation is used for 250 maintenance and growth respiration costs.

251
As in [42], we assumed the respiration to be proportional to the sucrose content. The former has a frequency of r m = 0.79 1 h . The frequency of sucrose loaded into phloem for growth has a frequency of η = 1.98 1 h . However, as observed in [45], the nightly starvation reduces the growth rate. To take into account this negative feedback, we reduced r g under the parameter γ that in the model describes the nightly efficiency of starch. Hence: The more γ increases, the more r g decreases. Over maintenance and growth respiration, sucrose is consumed for the transport of carbon along the plant and for the uptake of nutrients. [20] estimates costs of uptaking phosphorus as: with a sugar transport cost of 0.035 [−] (carbon consumed for each gram of carbon allocated in phloem). The cost of uptaking nitrogen (n c ) is assumed, in our model, to be the same of phosphorus (p c ), being the uptake mechanism similar among all chemical elements. Hence: Named u p [ µmolP gF W ·h ] and u n [ µmolN gF W ·h ] the phosphorus and nitrogen rates of uptake respectively, the whole uptake (r u m ) and transport (r t m ) costs are:

Further environmental contributes 252
Additionally, we should consider how deficiency or excess of nutrients in soil affect plant behaviour. For instance, [69] shows how in nutrient-poor soils, the secondary metabolism and the plant defence overcome the growth stimulus. For this reason, it is important to quantify nutrient excess or deficiency with respect to an optimal threshold. In [70] for free N medium the value of 12 µmolN cm 3 is used, and this value is used as optimal threshold in our model for nitrogen. [66] defined treatments in low P content up to 0.1 µmolP cm 3 , while treatments in high P content start from 0.2 µmolP cm 3 . Hence, we fixed the optimal value for phosphorus to 0.15 µmolP cm 3 . Therefore, the excess or deficiency will be measured with respect to the distance from the defined optimal thresholds and the amount of relative roots that are experiencing it: The increased defence-signal can be simulated limiting the conversion parameter λ sb , 253 and the increased metabolism will lead to greater nutrient stocks. We can simulate it 254 increasing n min and p min . Both N and P are primary for plant growth and a negative 255 feedback should be present to consider their combined effect.

256
The function R − (defined in Eq (38)) is used to reduce λ sb , regulating this way toxicity conditions, and it is symmetric with respect to n s and p s , with its maximum in (n s , p s ) = (1, 1). Moreover it is null if both nitrogen and phosphorus are missing in soil ((n s , p s ) = (0, 0)), while all negative values are also cut to 0. Hence: The function R + is instead used to increase nutrient demand, to regulate deficiency conditions in soil, and it is again symmetric and positive. It has one minimum R + (1, 1) = 1 and an arbitrary value was fixed for R + (0, 0) = 139. Hence: and n min = n min R + (n s , p s ) p min = p min R + (n s , p s ) where R + (n s , p s ) = 1 if min(n s , p s ) ≥ 1. possible negative effects) [71]. Moreover, the high nitrogen content in plant can increase 262 the secondary metabolism instead of growth [69]. Other examples come for the 263 phosphorus. According to [71], an high phosphorus content in soil can promote the 264 luxury uptake of P. The plant absorbs more P than it needs modifying P:Fe and P:Zn 265 stoichiometry ratios in plant. Consequences are the same of a Fe and Zn deficiency.

266
As case study, we wanted to include in our model the plant behaviour in case of luxury levels of P. To this aim, we should insert new equations for Fe and Zn dynamics and consequently more parameters to estimate optimal stoichiometry ratio. However, these data are not easily available in the literature. Hence, we introduced in the model a function simulating Fe deficiency when P content overcomes an optimal threshold. Iron is fundamental in many metabolic pathways and photosynthetic organs. In fact, [72] shows how Fe deficiency affects the development of chlorophyll, limiting the rate of photosynthesis. In the paper, the authors conduct some experiments pointing out the strength of this inhibition. It can induce up to a 50% decrement in chlorophyll accumulation. Moreover, [72] explains that is not possible to induce a greater decrement. Hence, if p s ≥ 1, the maximum photosynthesis is reduced depending on the amount of phosphorus excess, according to the following: If there is no excess, the photosynthesis does not experience any limitation 267 (p max h = 12.7). The more excess, the more reduction occurs till the 50%.
−a n ( n n + n min + n c u n n c u n + p h (1 − γ) + τ as ) In order to estimate them, we used the data available in [46]. In this work, the authors 286 grew Arabidopsis under 5 different photoperiods (4h, 6h, 8h, 12h, 18h of light) for 29 287 days assuming not limiting medium growth. Moreover, they measured sucrose and 288 starch content and the aboveground rate of growth.

289
Within our model, we assumed only one soil zone (b c r = +∞ and n s2 = p s2 = 0 µmol cm 3 ). 290 Moreover, as already said, no competition effects were taken into account

292
Of course, in estimating parameters we have to consider a strong dependence of natural processes by temperature. Since the influence of temperature was not explicitly modelled, we estimated different parameters for each photoperiod from [46], using a spline interpolation to approximate different day-lengths. For sake of simplicity, we tried to keep fixed the most of parameters. The best fitting was obtained assuming as variables only λ c , λ sni , λ sb , obtaining the following values:  We compared the results obtained by our model with the results extracted from the 301 biological data available in [46]. 302 We looked at the starch dynamic over different lengths of the day: 4, 6, 8, 12 and 18 303 hours (Fig 2A-C and Fig 3A-B). In each plot, red lines are piece-wise linear functions 304 interpolating data from [46], while green lines are the numerical results obtained from 305 our model. Blue dots represent the dawn (beginning of the day) and red ones the dusk 306 (beginning of the night). The starch grows during the period of sunlight, and decrease 307 during night with a good match between numerical and biological results. 308 Fig 2. Validation test in short days. Starch (from A to C) and sucrose (from D to F) dynamics for (A, D) 4h day length, (B, E) 6h day length, (C, F) 8h day length. The red lines represent the data obtained from [46], green lines are obtained from our numerical simulations. Blue dots are the beginning of the day and red dots are the beginning of the night. can be considered in average in agreement with biological data, even though less match 310 is found between peaks of dawn and dusk. In fact, the oscillations in sucrose content are 311 smoothed in the model. For instance, we can notice in 4h photoperiod (Fig 2D), instead 312 of having two sucrose peaks during night, the model is able to simulate only one higher 313 peak. A smoother behaviour is instead more evident in 18h of light ( Fig 3D); a closer 314 view on the oscillations for this photoperiod is reported in Fig 4A. Also in these plots, 315 red lines fit experimental data, green lines represent simulations, while blue and red 316 points are dawn and dusk respectively. Another defect concern the delay in the 317 numerical results with respect to the experimental data. For instance in the 8h 318 photoperiods (Fig 4B) there is a delay of 4h in the highest peak while in in the 18 319 photoperiods the delay is about of 6h. Anyway, the averaged dynamics well match 320 simulation and biological results. The reason for such delay can be imputed to the τ as function, which formulation 322 require more specific investigations (see file S1 in supporting information). Both the 323 delay and the smoother dynamic effects could imply an on/off control on τ as , with 324 respect to some thresholds. Moreover, by assuming all parameters depending on 325 photoperiods, it would be possible to further reduce the differences between model and 326 biological data. However, we considered this validation sufficiently sharp (see next 327 section 3.2), demonstrating the consistency of the model with respect to the literature. 328 Finally, we compared the relative growth rates (RGR). In [46], the RGR of fresh biomass of leaves is measured at the end of the night after 29 days. The mean values extracted from the literature are reported in Table 1 and compared with the RGR simulated by our model, obtaining an absolute error of 3.6283 · 10 −4 . To compute the RGR, we used the exponential formulation [73]: where t 0 is the initial starting time. To evaluate the robustness of our model, we compared the numerical results with data 333 extracted from several other papers. For instance, we tested the adaptation of the 334 model to different light periods. In [42], it is shown how plant changes its metabolism 335 and its sucrose content when day length changes from 8h to 16h (Fig 5).  [42]. Dynamic of sucrose content in sinks in Arabidopsis when passing from 8h to 16h of light (in top left plot) and vice versa (in the top right plot). The black line refers to wild type plants. 336 Since, we are not explicitly considering sucrose content in sinks, we simulated only 337 sucrose content dynamics passing from 8h light to 16h. The quantitative behaviour can 338 not be the same, but the qualitative behaviour is supposed to be consistent. In fact, our 339 model correctly predicts this dynamic as shown in   Just like in [42], in the day when photoperiod changes, plant adapts its behaviour in 341 an intermediate way, to fix it during successive days.

342
To verify the behaviour of γ we again referred to [42], in which γ varies along 343 different photoperiods and it is estimated according to Table 2.  Another important validation is about the N:P ratio. In [23] it is shown how plant 353 tends fast to optimal ratio when a missing nutrient is inserted in the environment where 354 the roots grow. We fixed the value 10 according to the same reference and a light period 355 Table 3. Numerical results for γ at different photoperiods. 6h 8h 9h 10h 11h 12h γ 0.7 0.63 0.6 0.56 0.52 0.49 of 12h. To test this behaviour, we fixed soil conditions and we varied the initial nutrient 356 ratio. Starting from a very small and an higher initial ratio, the behaviour tends to 357 approach the optimal ratio (Fig 7).  Interestingly, [76] pointed out the role of N:P ratio in quantifying nutrient limitation 359 on growth. Indeed, when we simulate different N and P soil contents, the plant reaches 360 different stoichiometry ratios, as showed in Fig 8 (the point highlighted refers to N:P 361 ratio in optimal nutrient soil content).

Fig 8. N:P dependence on soil nutrient contents
Moreover, when modelling the nutrient uptake, we assumed a dependence on the internal status instead of soil content. In [63], different values for both I max p and k p were proposed for different plants of Arabidopsis. We averaged among them in order to estimate the hourly phosphorus content: in poor p soils (p s = 0.0025 µmolP cm 3 ) and rich p soils (p s = 0.5 µmolP cm 3 ). Then we run our model and computed hourly uptake: In Table 4 there is the comparison between numerical results and our model, for 16h of 363 light and 7 days of growth, showing a maximum error of 3.52%. Finally, we looked at the effects of excess and deficiency of P on growth. In [77] are 365 reported two experiments on Arabidopsis, growing plants in 16h of light for 19 days 366 under two different P treatments. The first treatment consists of a limited P medium 367 (p L s = 0.006 µmolP cm 3 ), and the second consists of a saturated P medium (p H s = 1 µmolP cm 3 ).

368
The authors measured the total dry biomass (DW) and the root:shoot ratio each 4 days. 369 We used data from the 7 th day as initial point and we compared both fresh biomass and 370 root fraction at day 11 th , 15 th and 19 th . Since our model is parametrized on fresh 371 biomass, from [46] we estimated the relation DW= 0.08F W . Results, reported in Table 372 5, show an average error over the three periods of 0.028 ± 0.03 gDW in the rich soil and 373 0.00085 ± 0.00091 gDW in the poor soil, for the total biomass, and 0.43 ± 0.35 (%) in 374 the rich soil and 1.71 ± 1.36 (%) in the poor soil, for the root fraction. The errors in 375 biomass estimation, especially for the total biomass in rich soil, could be induced by the 376 approximation adopted for FW to DW conversion. In fact, the water content could vary 377 along the plant lifespan, affecting the relationship between dry and fresh weight. The evolution in time of the biomass and the root:shoot biomass ratio obtained from 379 our numerical simulation is also reported in

381
Up to now, soil conditions were stationary in time and nutrients were supposed to be only in the topsoil zone. In this section, we arbitrary fix b c r = 4 · 10 −5 gF W in order to study root system adaptation in heterogeneous soils. Results presented in the following refer to simulations of 20 days with 12h of light. We at first simulate root growth setting the topsoil rich in p, and the subsoil rich in n: From Fig 10A we can observe that no n is experienced in the subsoil zone (and 382 consequently no root is allocated) until the critical mass is reached. Then, plant starts 383 allocating resources to elongate in the subsoil zone. Of course, at the beginning the 384 more P reduces the N:P ratio that is immediately increased when the subsoil is reached 385 and the more nitrogen is available (Fig 10B). A more realistic simulation could be done assuming both nutrients initially placed in topsoil. But assuming that, day by day, nitrogen is transported by water from the topsoil to the subsoil: p s1 = 0.15 µmolP cm 3 , n s1 = 12 · 0.9 i−1 µmolN cm 3 , where the exponent i is the time (in days). Results are shown in Fig 11A. The less n s1 387 transported by water, the more n s2 accumulated in the subsoil.  To simulate a deeper nitrogen leak and then an higher topsoil zone missing n, the critical biomass b c r can be increased. Moreover, it can happen the death of the plant, if nitrogen is lost faster than the time needed by roots for soil exploration. This condition can be obtained by the absence of growth. Let's assume: In this case, the root biomass is not able to grow up enough to reach the subsoil, 389 although the greater investment in roots (different slopes in blue and green curves in Fig 390  11B). This way, the nitrogen is quickly exhausted and the plant is not able to survive. 391 We also verified the adaptive behaviour of plant in rich-soil zones (as reviewed in [55]): The plant used to grow in a more and more limiting topsoil, until the subsoil is reached. 392 The more nutrient availability stimulates the growth in the deeper soil zone (Fig 12A).

399
The adaptation of uptake signals a n (magenta line for nitrogen) and a p (light blue 400 line for phosphorus) is instead presented in Fig 12D. In the end, the balance between 401 uptake gain (u n and u p ) and nitrogen and phosphorus costs (C n and C p ) are shown in 402

403
All together, these plots allow us to appreciate a deeper simulation of plant 404 adaptation. In fact, at the beginning, the plant experiences both low n and p in the 405 topsoil. The limiting resources promote a greater uptake and a greater root priority. In 406 Fig 12C the root fraction quickly increases. Moreover, the phosphorus is more limiting 407 than nitrogen. In fact, in Fig 13 the averaged nitrogen balance is higher than the 408 averaged phosphorus balance (the difference is evident at the beginning of the graph). 409 This leads to a greater N:P ratio ( Fig 12B) and a smaller a n .

410
However, according to [55], differences in uptake costs and root growth costs define 411 the nutrient affinity improvement of plant or the increasing of root fraction. In the 412 phosphorus limiting case, being the soil exploration useless due to nutrient deficiency 413 (until the subsoil is reached), it is more convenient to increase phosphorus affinity a p 414 more than root fraction (light blue curve in Fig 12D and magenta curve in Fig 12C). On 415 the other hand, when the rich soil is reached, plant cannot immediately change its 416 strategies due to the concept of memory [18]. In fact, plots show no changes in 417 dynamics after having passed the topsoil boundary. Anyway, more phosphorus is now 418 available and the unmodified strategy leads to luxury uptake of phosphorus and a 419 critical reduced N:P ratio (Fig 12B). To stabilise the stoichiometry ratio, plant starts 420 reducing phosphorus affinity. Moreover, in rich soil zones, the uptake is more convenient 421 than root proliferation, so that root priority is less strong than a n . In the end, when 422 both the optimal ratio is reached and the memory is coordinated with the new soil 423 nutrient content, plant can stabilise its signals. biomass is reached after t = 343h. However, even after the critical biomass, plant does 429 not stop growing (Fig 14)but the sucrose is not enough to sustain the actual growth 430 (Fig 15A). This is the signal that later will reduce the growth. Both leaves and roots 431 decrease with respect to a plant not experiencing the competition (Fig 15B and 15C).

432
Of course, there are no reasons why affinity should change (for example in Fig 15E we 433 report the phosphorus affinity a p being similar also the a n dynamics). Letting 434 conditions unchanged and simulating 50 days, a new equilibrium is reached by the plant 435 ( Fig 15E).  According to [78] and [79], the competition among adjacent roots occurs when the 438 respective rhizosphere overlaps the other one, reducing the uptake for unit of root mass. 439 This is consistent with our modelling approach (see section 2.6.1). For this evaluation, 440 we assumed only one soil zone, b max r1 = 0.3gF W and, initially, no leaf competition. We 441 performed simulations for 20 days at 12h of light, and the results are shown in Fig 16. 442 As expected, plant promotes leaf growth instead of the root one.  In the former case (greater leaf critical biomass), the strategy is dangerous because 447 favouring leaves than roots, roots can not absorb enough nutrients to sustain the faster 448 leaf growth. At the beginning, the total biomass is greater with respect to the no 449 competition framework (red curve in Fig 17A).  After few days (25 instead of 20 days) the plant dies (red curve in Fig 17B). In the 451 second case, the leaf competition happens when the plant can still sustain its 452 development. The growth still hold but in a slower way and the plant does not die (red 453 curve in Fig 17C). experiences the competition inducing an unbalance between nutrients: the N:P ratio is 457 not anymore optimal (Fig 18A). After having reached the subsoil, the plant promotes 458 the root growth in the subsoil (Fig 18B and 18C) and stabilises again the optimal ratio. such as photosynthesis, uptake or water transport; they are devoted to the description 468 of the chemical interaction among hormones; or limited to a single nutrient.

469
To have a glance on the internal communication driving the growth, and differently 470 from the state of the art, we proposed here a model that integrates at the same time several processes (photosynthesis, starch degradation, nutrients uptake and 472 management, biomass allocation, maintenance), different signals (circadian clock and 473 growth stimulus) and nutrients (sucrose, nitrogen and phosphorus).

474
The overall system extends results from previous models ( [20,41]) and proposes an 475 approach, that could be addressed to biologists, for analysing interactions not yet 476 investigated. For example, thanks to the model it is possible to evaluate how 477 photosynthesis intensity and sucrose production rate are reduced or increased according 478 to different stoichiometry ratio values, or which process is mainly, and earlier, affected 479 by nightly starvation. By evaluating the dynamics of the priority signal f r when the  Comparing our model with [20,42], a new relation has been defined for the growth, 489 which does not depend anymore just on nutrients availability, but it takes into account 490 also starvation and efficiency. Here, by efficiency of an organ (which could be a leaf or a 491 root), we mean the ability of that organ to provide the amount of resources needed for 492 the correct functioning and growth. For example, the roots efficiency is the ability to 493 provide both nitrogen and phosphorus such that nutrient levels are enough to sustain 494 hourly costs of nitrogen (C n ) and phosphorus (C p ). On the other hand, leaves efficiency 495 is the ability to produce enough sucrose and starch so that all respiration costs (r m , r u m , 496 r t m , r g ) are covered by the available sucrose. Hence, as reviewed in this paper, plant 497 tries to never reach minimum nutrients thresholds, while keeping monitored level of 498 stores to avoid consuming energy unnecessarily.

499
What the model can show, by analysing these dynamics, is a dense network of 500 signals and feedback that enforce the concept of intra-plant communication allowing a 501 deeper understanding of how this internal network is connected.

502
Nevertheless, the model can be used to better define biological concepts like 503 performance, nutrient efficiency and sink priority.

504
In [24], performance is defined as "the ability to acquire resources and survive in the 505 presence of competition or in stressful environments where physiological limits are 506 reached". The problem in measuring the performance refers to the problem of defining 507 the parameter to be checked in stressful conditions: either the biomass, sucrose 508 production, nutrient uptake or photosynthate allocation are important parameters that 509 can provide an idea about the ability of the plant to survive. Our model allows to easily 510 estimate the plant performance by simulating limiting scenarios and measuring sucrose, 511 nitrogen and phosphorus levels with respect to the critical thresholds.

512
In [80], the efficiency in using nutrients is generally defined as: the "measure of how 513 well plants use the available mineral nutrients. It can be defined as yield (biomass) per 514 unit input (fertiliser, nutrient content)". A similar definition can be found in [25] for 515 efficiency in the use of water, and related to photosynthesis efficiency. 516 Clearly, the biomass parameter is not enough to measure the efficient use of nutrient; 517 in fact, just considering this parameter, a simple imbalance in roots and leaves biomass 518 could imply a not efficient use of resources, whereas it is not always the case [81]. In our 519 model, the efficiency emerges from relationships among uptake and nutrient costs (from 520 Eq (29) and (30)). In fact, it is possible to note a nutrient efficiency in the first ratio, 521 relating uptake (u n and u p ) with costs (C n and C p ), and a photosynthesis efficiency in 522 the last ratio, which relates uptake and photosynthesis. In addition, Eq (13), relating 523 sucrose and its critical thresholds, gives an idea about sucrose efficiency.

524
Finally, [82] defines the strength of sinks, or sink priority, as "the competitive ability 525 of a sink organ to import photoassimilates, and depends on both its physical (size) and 526 physiological (activity) capabilities". While the size is a clear parameter and it can be 527 easily related to the biomass of the organ, defining the activity of a sink is more 528 difficult. In [83], the term sink activity is referred to all processes involved in using 529 resources, from the synthesis of new tissues to the maintenance of the already existing 530 ones. This activity can be quantified as the net assimilation of resources in tissues 531 minus losses due to respiration and exudation. However, in [84] it is noted how the 532 previous definition of sink strength is not enough to explain how sucrose is allocated in 533 the plant. In fact, processes affecting osmotic concentrations in sinks need to be taken 534 into account for resources allocation. Furthermore, in [83], is highlighted the importance 535 of associating the source-sink interaction with the whole plant growth, even if no hint 536 about how to do it is given.  swarm robotics [23], for resolution of combinatorial problems [85]).

553
The mechanisms of plant communication are being studied through self-recognition 554 cues, underground interactions with fungi or other plants, volatile organic compounds, 555 electric and chemical signals ( [86][87][88]). In this paper, we investigated in the intra-plant 556 communication by analysing the possible cues activating the adaptive growth responses 557 in a single plant and describing the dynamics of such internal processes and signals with 558 a mathematical model. Our assumptions and formulations have been deeply analysed in 559 comparison with other models from the state of the art, highlighting the adopted 560 extensions and improved potentialities.

561
Main difference with other models is having included several different elements (e.g. 562 photosynthesis, starch degradation, multiple nutrients uptake and management, biomass 563 allocation, maintenance) and their dynamics all at once considering their interactions 564 and effects on the growth, thus showing the potential intra-communication in plants; 565 whereas, strength of our model lies in having based assumptions and formulations on 566 biological evidences and laboratory experiments collected from the state of the art. This 567 approach allowed to evaluate several parameters from the papers thus reducing the 568 number of free variables, and at the same time, to maintain a high level of details with 569 good agreements with the biological model.

570
To test the robustness of our model, we compared the biomass, sucrose dynamics, 571 photosynthesis partitioning rate, uptake strength and sinks priorities not using custom 572 experiments, but encouraging the comparison with different independent published 573 biological data. Indeed, we validated the model by comparing experimental results with 574 results obtained in our simulations, reproducing conditions of growth similar to the 575 biological counterpart. All validation tests show high accuracy in results and very small 576 errors (e.g., we obtained 3.52% as maximum relative error in nutrient uptake 577 comparisons (Table 4) and 3.44% as maximum percentage error for biomass estimation 578 in low phosphorus soil tests (Table 5) ); where, the errors obtained are likely induced by 579 divergences of initial conditions between biological and numerical experiments. However, 580 both the qualitative behaviour and the quantitative values agree with data in literature. 581 Moreover, we simulated different growth scenarios. Interestingly, in stress conditions, 582 e.g. when the saturating thresholds were reached, the dynamics showed a fast 583 adaptation of the plant with a self-optimisation of internal resources and feedback 584 stimuli, without the need of a a-priori optimisation assumption. The same self-adaptive 585 behaviour was evident when branching vs elongation was analysed. This optimising 586 behaviour is realistic but difficult to be biologically demonstrated. However, differently 587 to some models [4,19] that directly start by assuming the existence of an internal 588 optimisation function which drives plant dynamics, in here we demonstrated it to be a 589 consequently behaviour emerging from the dynamics of internal processes. Therefore,  [23]). For instance, the equations, identically replicated on each robot, can 597 describe the behaviour of the robot and represent the ability of robotic roots to colonise 598 a certain area (plant branching) or to growth further in exploration (plant elongation), 599 according to the internal signals (into the colony of robots) and the measurements from 600 the environment.

601
Future steps will include from one side, the implementation of the proposed 602 behaviour in robotic solutions, while from a biological perspective several improvements 603 can be adopted. For instance, the saturating thresholds can be replaced by death terms 604 that need to be estimated via laboratory experiments. Stem and fruit dynamics can be 605 introduced for a more complete description. Moreover, the concepts for resources 606 allocation, analysed in this paper only in roots, could be used to study resources 607 partitioning in leaves. Additionally, from a mathematical point of view, partial 608 differential equations can replace root dynamics to analyse spatial diffusion of roots.
609 Moreover, the model is easily embeddable with more detailed formulations for 610 photosynthesis, light, temperature, starch degradation, respiration or exudation in order 611 to make the model more and more detailed.