A mathematical model explores the contributions of bending and stretching forces to shoot gravitropism in Arabidopsis

Abstract Plant shoot gravitropism is a complex phenomenon resulting from gravity sensing, curvature sensing (proprioception), the ability to uphold self-weight and growth. Although recent data analysis and modelling have revealed the detailed morphology of shoot bending, the relative contribution of bending force (derived from the gravi-proprioceptive response) and stretching force (derived from shoot axial growth) behind gravitropism remains poorly understood. To address this gap, we combined morphological data with a theoretical model to analyze shoot bending in wild-type and lazy1-like 1 mutant Arabidopsis thaliana. Using data from actual bending events, we searched for model parameters that minimized discrepancies between the data and mathematical model. The resulting model suggests that both the bending force and the stretching force differ significantly between the wild type and mutant. We discuss the implications of the mechanical forces associated with differential cell growth and present a plausible mechanical explanation of shoot gravitropism.

Pioneering data studies of gravitropism were based on the development of image analysis tools to quantify the tip angles of shoots and roots (Coutand et al., 2007;Kutschera, 2001;Philippar et al., 1999). However, such a simple angle extraction o en fails to account for the actual morphology of bending organs (especially for shoots) because the shoot has a continuous filament-like structure and bending involves changes in curvature and length that are distributed along the entire shoot [the limitations of quantifying only tip angles are summarized in Moulia and Fournier (2009). To solve this problem, a new technique that quantifies continuous displacement by adding a distinctive randomized pattern to an object was developed to track the growth of roots [KineRoot, Basu et al. (2007)]. Recently, KymoRod was developed to recognize root angle, curvature and relative growth rate without patterns (Bastien et al., 2016). e data approach is now widely applicable with high precision but omits mechanical information, such as stretching and bending force, that can only be extracted using model approaches.
Researchers have constructed a variety of mathematical models to describe the shoot's morphology and bending during gravitropism (Bastien et al., 2013;. Early models were based on active processes where the curvature changes depending on the horizontal angle of the stem segment relative to the ground (gravity sensing) and the curvature of the stem segment (straightening or proprioception) (Bastien et al., 2014). Although some of these models have analytical solutions, analytical studies are not applicable to all types of shapes a er bending. Recently, based on the active process mentioned above (Bastien et al., 2014), a numerical method that ensures the minimization of elastic energy has been developed including the effects of elasticity; this method explains how stationary shapes are affected by the inclination angle and curvature of the shoots (Chelakkot & Mahadevan, 2017;Moulia et al., 2019). One of these studies (Chelakkot & Mahadevan, 2017) revisits the concept of 'morphospace' , a well-known methodology in evolutionary and anatomical studies (Eberle et al., 2014;Polly & Motz, 2016) whereby representing a possible form, shape and structure of an organism only by a few characteristic parameters (e.g., width, height, growth rate, and so forth). is type of modelling approach is important because it incorporates the mechanical aspects of shoot gravitropism.
Despite these advances in data and model approaches, the role of mechanical forces during bending of plant organs is poorly understood. e purpose of this study was to elucidate the mechanics of gravitropism using mathematical model based on actual data. First, we analyzed the morphological differences between wild-type Arabidopsis and the lzy1 mutant during shoot bending. Next, we combined the data with the previously reported model in Chelakkot and Mahadevan (2017) to extract mechanical forces applied to the shoot during gravitropism. Finally, we formulated the biological implications of these mechanical forces and discussed possible experiments to verify them.

Maximum curvature a er 140 min of bending is significantly lower for the lzy1 mutant than for the wild type
We quantified the morphological state of shoots before and during bending, by reanalyzing gravitropic responses which had been reported previously (Taniguchi et al., 2017). Time-lapse images of wild type and lzy1 mutant inflorescence stems bending in response to gravity are shown in Figure 1a,b. We used ImageJ to extract the centerline of each shoot from the original images and constructed a continuous description of the centerline by spline interpolation (see Section 4). Based on the continuous curve, the evolution of shoot lengths and curvatures were calculated for 16 wild-type and for 15 lzy1 mutants. ree representatives of each genotype are shown in Figure 1c,d (see all examples in Figure SI2). We set the origin of the curvilinear abscissa (mm) at the base of the shoot, with values increasing with distance along the centerline from the base to the tip. e curvature at the middle of the wild-type shoot changed to positive values (red) a er 40-50 min ( Figure 1c) but was lower for lzy1 mutant shoots (light-red) (Figure 1d).
For a morphological description, we produced a morphospace ( Figure 1e) using the extension ratio of shoot length (corresponding to shoot strain) and the maximum curvature at 140 min as characteristics of shoot shape. e extension ratio was calculated by comparing the initial and final shoot lengths (L t e − L t s ) L t s , where L t s (L t e ) denotes the shoot length at the starting (ending) time t s (t e ) a er placing the shoot in an almost horizontal position. We note that the morphospace can be described using other characteristics (e.g., inclination angle, timing of bending, and so forth), but we chose the two characteristics above for the following reasons. e extension ratio enables us to check the consistency of our results with the reported observation that the lzy1 lzy2 lzy3 triple mutant does not show significant growth impairment (Taniguchi et al., 2017, which we assumed would also be the case for the lzy1 mutant. Using the maximum curvature, it is possible to check if there is a significant difference in the degree of bending between the wild type and the lzy1 mutant, as reported in Taniguchi et al. (2017). Consistent with these previous results, the morphospace indicates that the extension ratio of the lzy1 mutant is not significantly different from that of the wild type, but the maximum curvature at 140 min is significantly lower in the mutant than in the wild type ( Figure 1f).

A mathematical model captures typical shoot bending events in the wild type and the lzy1 mutant
To understand the mechanics underlying the different morphologies of the wild type and the lzy1 mutant, we reimplemented the mathematical model proposed in Ref. Bastien et al. (2013Bastien et al. ( , 2014, Chelakkot and Mahadevan (2017) (see Section 4). We note that the mathematical model itself in this study is exactly same as the model in Chelakkot and Mahadevan (2017) but the fitting of the experimentally obtained bending dynamics data to this model is novel in the context of gravitropism. In this model, the plant shoot is modeled as a growing elastic rod under gravity. e rod is separated into growing and nongrowing zones. In the growing zone, the intrinsic (or spontaneous) curvature κ * and the natural length of the stem segments change as functions of time during gravitropic and proprioceptive responses . At each time step, the shoot configuration is determined through force and momentum balance equations of the rod under the target geometry. By contrast, in the nongrowing zone, we do not change the target geometry. To compute the shoot shape of the model, we discretized the shoot centerline into a set of particles (vertices) connected by elastic springs (Figure 2a le ). e positions of the particles are determined by the balance of stretching, bending and gravitational forces (Figure 2a right). We repeated the updates of the target geometry and the current configuration to compute the shoot shapes.
is model system has two characteristic dimensionless parameters: the growth-sensitivity parameter S and growth-elasto-gravity parameter ε. e parameter S represents the ratio of the size of Morphospace diagram with maximum curvature and extension ratio for data and model for the wild type (blue) and for the lzy1 mutant (green). The square plots represent the averaged values for the wild type (blue) and for the lzy1 mutant (green). Error bars are standard deviations. (f) Boxplot of the maximum curvature at 140 (min) showing a statistically significant difference in maximum curvature between the wild type and lzy1 mutant with t-test p < .005 (le panel), and not in the extension ratio (right panel). the growing zone l g relative to the sensitivity length l s , that is, S = l g l s . e parameter ε is the ratio of the elasto-gravity length l e = (B ρg) 1 3 and the size of the growing zone, that is, ε = l g l e (see also Section 4).
Based on the morphospace in Figure 1e, we systematically searched through values for the dimensionless parameters S and ε to find those that resulted in the minimum deviation between the model and data. We note that the search was done with fixed l g with the experimentally determined initial length, and with various l e for ε by changing mass density per unit length, ρ and with various l s for S by changing both gravi-sensitivity β and the proprioceptive sensitivity γ (see Section 4). To extract the typical bending behaviour of the wild type and the lzy1 mutant, we identified the best-fitted parameters with the averaged values for the wild type (wt-ave model) and lzy1 mutant (lzy1-ave model). e extracted values of gravitropic sensitivity β in wild type (β = 2.14 ± 1.95) was higher than those in lzy1 mutant (β = 1.38 ± 1.30), reflecting the defects in DCG in lzy1 mutant. On the other hand, the extracted values of proprioceptive sensitivity γ in wild type (γ = 0.08 ± 0.07) was as same level as those in lzy1 mutant (β = 0.11 ± 0.08), indicating the same level of proprioceptive strength.
e averaged models exhibit similar bending behaviours as the real experimental data, where the wild type bends to a greater extent than the mutant does, as shown in Figure 2c (wild type) and Figure 2d (lzy1 mutant). While the model does not reproduce the shapes of the shoot perfectly due to the variation of the initial conditions in data (see Figure SI5), the spatio-temporal maps of the curvature were not significantly different from the data. e model predicts a spatio-temporal trend in curvature increasing from the tip to the bottom of the shoot overtime, and those decreasing at tip a er 100 min (Figure 2e,f), which are consistent with the real curvature data (Figure 1c,d).

The shoot models provide information about the mechanical forces involved in bending
e configuration of the model shoot is realized by the balance between the forces derived from the biological active process (the stretching and bending forces), andexternal gravitational force. We note that our model shoot is inextensible and unshearable filament as the model in Chelakkot and Mahadevan (2017) or other models of growing elastic rod in Goriely (2017) at every mechanical equilibrium configuration. e stretching and bending forces we discuss in the next section originate from changes in natural length and intrinsic curvature, respectively. To clarify the role of the stretching and bending forces, we decomposed the total displacement of the shoot → u at each time step into the individual displacements from the stretching force → u s , the bending force → u b and gravity → u g in the wt-ave model (Figure 3a,b) and in the lzy1-ave model (Figure 3c,d). We note that the displacement from gravity was negligible both in wt-ave and lzy1-ave model.
In both the wt-ave and the lzy1-ave models, the shoot tip deforms greatly, but the displacement at the base is small (Figure 3a1-d1). We noticed that the stretching and bending forces work differently at the tip and at the base. e direction of the stretching force near the tip is opposite to the bending direction, implying that the stretching force near the tip suppresses bending (a proprioceptive effect), as shown by the blue-filled arrow in Figure 3b2. However, bending near the tip is promoted by the direction of the bending force, as shown by the red-filled arrows in Figure 3b3. Near the base, by contrast, the stretching force promotes bending, whereas the bending force suppresses bending (red and blue open arrows in Figure 3b2,b3, respectively). Because there is little morphological change at the base, it is expected that the effects of these two forces on bending are comparable in magnitude, resulting in no significant change at the base.
e lzy1-ave model exhibits spatio-temporally similar behaviour, although the colour diagram shows less displacement and smaller changes in stretching and bending forces compared to the wt-ave model (Figure 3c,d). Accordingly, the maximum stretching force and maximum bending force are significantly weaker in the lzy1-ave model compared to those in the wt-ave model (Figure 3e,g). Quantitatively, the magnitudes of these forces are about 50-60% smaller than those in the wt-ave model. at means that the defects in the signal transduction from the gravi-sensing in lzy1 mutant are reflected in the weaker magnitudes of stretching and bending force.

Biological implications for gravitropism, proprioception and shoot axial growth
We have clarified that the bending and stretching forces have different effects on shoot morphology during bending. e remaining question is whether the decomposed bending and stretching forces can be translated into biological counterparts. e biological process underlying the change of the curvature is differential cell growth (DCG), which is realized by a combination of the gravi-proprioceptive response and the shoot axial growth as discussed below. e bending force controls the degree of gravity-free bending and is associated with both gravitropism and proprioception through the change in the intrinsic curvature κ * . Gravitropism promotes bending, whereas proprioception suppresses bending, as illustrated in Figure 4a. e former and latter depend on the inclination angle and the curvature, respectively. We call this combined effect gravi-proprioceptive response (GPR). e gravitropic response becomes significant when a stem segment has a low inclination and low curvature because, in this case, gravitropism is stronger than proprioception, resulting in the promotion of bending (Figure 4a top). By contrast, a high inclination and high curvature suppress bending because proprioception is dominant (Figure 4a bottom). erefore, the bending forces for the wild type are significant at the start of the bending process and become smaller, starting from the tip, as the shoot bends (Figure 3a3,c3). e stretching force is determined by shoot axial growth and is governed by both the growth and the curvature of the stem segments. As the shoot grows, it changes the intrinsic curvature κ * and the current curvature κ is governed by mechanical balance. We call this combined effect growth-curvature-induced bending (GCB), which is suppressed as the shoot grows because the corresponding tensional forces act in the direction opposite to the bending direction, thereby suppressing bending (Figure 4b top). By contrast, GCB is enhanced in the absence of growth because the contracting forces act in the same direction as bending (Figure 4b  bottom). For the wild type, therefore, the stretching force is less significant at the start of bending, and is reinforced around the base, as shown in Figure 3a2,c2.
As a consequence of the two bending effects, the dynamic shoot morphology (i.e., curvature) is determined. Assuming the shoot is a solid material with thickness 2δ and relative elongation growth ratesε outer andε inner for the outer flank and inner flank, respectively (Figure 4c), the degree of the DCG is defined as ∆ (s,t) = (ε inner −ε outer ) (ε inner +ε outer ) [see Ref. Bastien et al. (2014)].
We then quantified the DCG asL 0 ∆ (s,t) =ε inner −ε outer . DCG is negative if the outer cells grow faster than the inner ones; otherwise DCG is positive. e DCGs for both the wild type and lzy1 mutant are negative when the shoot bends upward, but a er 100 min, when the shoot straightens, DCG becomes positive near the tip (Figure 4d). ese behaviours of the DCGs are consistent with the fact that displacement becomes large at the tip early in the bending process, then diminishes a er 100 min (Figure 3a1,c1). erefore, the mechanics of the stretching and bending force discussed above could be interpreted as the GPR and GCB, respectively, and the consequence of the two effects can be observed in DCG.

Future perspective: prediction for growth at the cellular level
e mathematical model provides a quantitative assessment of the differential growth at the cellular level. At the initial bending stage, the relative elongation growth rate of the outer flank at tip in the wild type wasε outer = 2 × 10 −4 min −1 and that of the inner oneε inner = −2 × 10 −4 min −1 (see details forε inner andε outer in Figure SI6), implying that the inner flank shrinks as the shoot grows. is counter-intuitive prediction at the cellular level may be experimentally clarified by measuring the growth rates of cells at the inner and outer franks for wild type and for several mutants with defect in building secondary cell wall (Kubo et al., 2005;Mitsuda et al., 2008) and with defect in posture control of shoot (Okamoto et al., 2015).

Discussion
In this study, we extracted the detailed morphologies of wild type and lzy1 mutant Arabidopsis shoots during bending. We confirmed that their maximum curvatures were significantly different, whereas their extension ratios were not. Subsequently, we fitted our mathematical model to the bending events and searched the best fitted dimensionless parameters, S and ε. We then constructed mathematical models corresponding to the averaged bending behaviours of the wild type and the lzy1 mutant and confirmed that the reconstructed bending events in the models were similar to the actual bending events in the data. With this setup, we could explore the mechanical information associated with stretching, bending and gravity, which is not observable in experiments.
e mechanics underlying the shoot bending events can be summarized as follows. e biological process underlying the change of the curvature is DCG, which can be interpreted as a combined result of the GPR and the GCB, which correspond to the bending force and the stretching force, respectively. e bending force (GPR) becomes larger if gravitropism is more significant than proprioception. Shoot growth stretches the shoot axially and induces another bending deformation, called GCB. As a consequence of the two bending effects, we can understand the dynamic shoot morphology in terms of mechanics. ese mechanical aspects of bending are the underlying mechanism hidden in the experimental description of shoot gravitropism.
For the improvement of the mathematical model, more realistic growth or bending models will need to be developed (e.g., the effect of lignification of secondary cell walls (Chelakkot & Mahadevan, 2017) and the temporal delay in gravitropic response (Agostinelli et al., 2020;Chauvet et al., 2019). For the improvement of fitting between data and model, in addition to our restricted two parameters (maximum curvature a er 140 minutes and extension ratio), the spatio-temporal curvature will need to be made correspondence for further understanding. is type of interdisciplinary study highlights the potential of connecting experimental observations with theoretical models and leads to a richer understanding of shoot gravitropism.

Data of gravitropism
We reanalyzed the data of gravitropic responses of wild type and lzy1 mutant which had been reported in our previous study (Taniguchi et al., 2017). We used A. thaliana Columbia-0 as the wild type and lzy1 (GABI_591A12). Surface-sterilized seeds were sown on MS plates [1 × Murashige Skoog salts, 1% (w/v) sucrose, 0.01% (w/v) myoinositol, 0.05% (w/v) MES (2-(N-morpholino) ethanesulfonic acid) and 0.5% (w/v) gellan gum; pH 5.8], incubated in darkness at 4 ○ C for 2-3 days, grown at 23 ○ C in a growth chamber under continuous light for 10-14 days, transplanted to soil, and grown under continuous light. Plants with primary inflorescence stems of 4-8 cm in length were gravistimulated by placing horizontally under nondirectional dim light at 23 ○ C a er 1 hr of preincubation. Photographs were taken at indicated times. 15 or 16 individuals of wild type or lzy1 mutant were tested, respectively.

Extraction of angle and curvature of shoot in data
To extract continuous angle information of the shoot, we used ImageJ to detect the positions of the shoot at intervals of approximately 0.2 mm. Subsequently, we interpolated the detected points (approximately 20 points) by a third-degree spline interpolation ( Figure SI1).

Mathematical model
We review the mathematical model proposed in Chelakkot and Mahadevan (2017), where gravitropic kinematics and mechanics are combined. e centerline of the shoot is modelled by a curve, whose position vector at time t is given by r (s,t) = (x (s,t),y (s,t)) . Here, s represents the arc-length of the shoot, satisfying 0 ≤ s ≤ L(t), with the total length L(t) at time t. We introduce the bending angle θ (s,t), as the angle between the local tangent and y-axis. Here, the tangent vector of the shoot is described ast = (sin θ, cos θ). We clamp the basal end (s = 0) as θ (0,t) = −π 2,x (0,t) = y (0,t) = 0, and the apical end (s = L) is set to be free for the force and moment at any time t. e shoot grows uniformly at speedL 0 in the growing zone of length ℓ g from the apical end. In the growing zone, the local intrinsic curvature of the shoot κ * (s,t) is modified by the active bending of the shoot, such as the local gravitropism and the proprioception for the current shoot configuration θ (s,t). e change of κ * in time is described as where β and γ are gravitropic and proprioceptive sensitivities, respectively (Bastien et al., 2014). e shape of the model shoot is determined by the moment and force balance of the rod with the intrinsic (inhomogeneous) curvature at every time t (i.e., the speed of the growth is assumed to be sufficiently slower than that of the mechanical relaxation). M (s,t) and F (s,t) = (H (s,t),V (s,t)) are, respectively, the internal moment (around z-axis) and force acting on the cross section of the shoot at the arc-length s and time t. e moment and force balance equations are given by ∂M ∂s with the constitutive law given by M = B κ − κ * . We have introduced the bending modulus B, the mass density per unit length ρ and the gravitational acceleration g. ree different length scales exist in this model system: the characteristic sensitivity length ℓ s = γ β, the length of the growing zone ℓ g and the elasto-gravity length (persistent length of elastic bending against gravity) ℓ e = (B ρg) 1 3 . us, the two critical dimensionless parameters are the growth-sensitivity parameter S = ℓ g ℓ s and the growth-elasto-gravity parameter ε = ℓ g ℓ e .
To simulate the shoot morphology, we discretized the centerline into a set of connected particles as r (s,t) → r i (t) (i = 1,2, . . . ,N), where i represents the index of the node (vertex). e position of the node r i (s,t) is updated based on the force and moment balance equation (2). e elastic force at the ith node is computed from the stretching of the bond b i = r i+1 − r i and the angle of adjacent bonds φ i . e stretching and bending potentials are respectively given by e indices d i and φ * i are the natural length of the bond and the natural angle of adjacent bonds (calculated from κ * ), respectively, which will be updated by the growth rule and equation (1). As assumed in Chelakkot and Mahadevan (2017), we chose the rod elasticity parameters such that the rod is practically inextensible. e typical dimensions of the shoot in our experiments are the initial shoot length L 0 = 40 mm (natural length of the bond d ∼ 1.33 mm) and the radius δ = 0.5 mm. e shoot has the Young's modulus of E y ∼ 10 MPa and the bending modulus = E y πδ 4 4 ∼ 5×10 −7 Nm 2 . e typical stretching force f * s and bending force f * b of an elastic rod that has the same E y and D are estimated as f * s = E y πδ 2 ∼ 8N and f * b = D L 2 0 ∼ 3 × 10 −4 N, respectively, from which we find f * s ≫ f * b . In the discretized simulation, this condition is realized by setting EL 2 0 B ≫ 1. We adapted the elasticity parameters E and B computed from those of a shoot via E = E y πδ 2 d and B = D d, which satisfies f * s ≫ f * b . In this paper, as the extension ratio of wild type and lzy1 mutant were similar (see Figure 1), we assumed that the elastic properties E and B for wild type are the same as those for lzy1 mutant. e total force applied on the vertex i is then given by where U g represents the gravitational potential. In the main text, we assign the stretching, bending and gravitational forces as (Chirico & Langowski, 1994;Sano et al., 2017). We note that the increase of bending rigidity with time associated with lignification of the secondary cell wall discussed in Chelakkot and Mahadevan (2017) is not considered in this study because it has a minor effect on the stationary morphology. e experimentally observed temporal delay of gravitational sensing discussed in Agostinelli et al. (2020), Chauvet et al. (2019) is also not considered because our focus is not the detailed fitting of the bending process but the extraction of mechanics with qualitatively consistent bending behaviour.
To clarify how parameters in the model affect the shoot morphologies, we changed the parameters in the model. rough an extensive parameter search, we showed how the following six parameters affected the shoot morphology ( Figure SI3). Note that though these parameters affect the shape morphology, the variety of shapes can be summarized by the two dimensionless parameters S and ε. e first three paramters are growth-related parameters l 0 , l g andL, defined as the initial length, growing zone measured from tip and the linear growth rate, respectively (illustrated in Figure  SI4). e other parameter l e is the elasto-gravity length arising from the competition between gravity and elasticity. e remaining two parameters, β and γ, are curvature-related parameters defined as sensing strength of inclination angle and curvature, respectively.

The exhaustive search for the best-fitted parameters S and ε
To detect the best-fitted parameters S and ε of the model (S model ,ε model ), we introduced a deviation D between data and model with fixed (S data ,ε data ) as follows To search the parameter ε model , the mass density per unit length ρ was systematically changed. To search the parameter S model , as the gravi-sensitivity β and the proprioceptive sensitivity γ are interrelated, we firstly searched both parameters β and γ and found a minimum deviation of D around the range S ∈ [2.0 ∼ 4.0]for all examples of the wild type and the mutant (see all the detected values (β model ,γ model ) and detected S in Figure SI7). As the parameter S is expected to be constant, and we detected that the parameters β and γ have a relationship γ model ∼ 0.078β model for all the data, we fixed the average S model ≃ 3.15 = <l g > 0.078 , < ⋅ > is sample average as an representative value and searched only the parameter β for each individual. e data that support the findings of this study are available from the corresponding author, Satoru Tsugawa, upon reasonable request. We share the simulation code that calculates the shoot morphologies on web browsers at github (https://satorutsugawa.github.io/linevertex.html).