Shapes of leaves with parallel venation. Modelling of the Epipactis sp. (Orchidaceae) leaves with the help of a system of coupled elastic beams

The model

The modelling of soft materials, in particular organic ones, is a notoriously difficult problem. Indeed, any functional organic tissue usually contains millions of elements, namely the cells and furthermore, such tissue is an arena of very complex chemical reactions combined with diffusion. Needless to say, the above statements also hold in the case of plant leaves. There are at least two possible theoretical approaches to the modelling of leaves. Firstly, one can think about simply simulating the shapes of the leaves and secondly, one can attempt to start (almost) from the first principles and include the cellular structure of the leaf tissue. With regard to plant morphology in general, including that of leaves, tremendous progress has been recently achieved within the framework of L-systems. Furthermore, the almost realistic pictures in two dimensions can be obtained by application iterated functions in a similar way to the construction of the Mandelbrot set (Bird & Hoyle, 1994). Bird and Hoyle’s simple approach has, however, nothing to do with any physical, chemical and biological properties of living organs. A more sophisticated approach has been developed using iterative functions by Runions et al. (2005) for the modelling of leaf venation. A powerful method to model biological pattern formation is the numerical solution of various reaction–diffusion equations (Koch & Meinhardt, 1994). A method to study general biological systems which exhibit various line-like and net-line structures using coupled reaction–diffusion equations was developed in the seventies (Meinhardt, 1976). Markus, Böhm & Schmick (1999) made an interesting attempt to model reaction–diffusion equations using cellular automata and then applying them to vessel morphogenesis including leaf venation. More recently, the shapes of leaves have been modelled by this method by Franks & Britton (2000). Another interesting contribution based on hydrodynamics rather than reactions and diffusion has involved the numerical solutions to the Navier–Stokes equations (Wang, Wan & Baranoski, 2004). Another starting point in leaf simulation has been worked out in Liang & Mahadevan (1999). The shapes of long leaves have been obtained using the standard theory of elasticity. Saddle-like midsurface and rippled edges of leaves have been interpreted as being caused by elastic relaxation through the bending that follows differential growth. The first-principles line of dealing with the problem of leaf modelling is still in statu nascendi; there is, however, promising work on the simulation of plant tissues (Ghysels et al., 2009; Van Liedekerke et al., 2011) based on micromechanical models of single cells (Van Liedekerke et al., 2010). Our approach lies in between a purely phenomenological approach and one based on the cell micromechanics. We realise that in the case of long leaves with parallel venation the anisotropy and inhomogeneity are essential, and models which assume isotropy and/or homogeneity must fail. This contrasts our model with that of Liang & Mahadevan (1999). However, we agree that elasticity theory is an appropriate framework for the description of leaf morphology. Indeed, regardless of what happens in the cells and what interaction between cells may be detected, the result must always involve changes in purely mechanical quantities which characterize a solid material. A continuum-mechanical model makes it possible to relate the elastic properties of the model with those of real leaves. We note here that the experimental research on the mechanical properties of leaves appears to have only just begun, and we are aware of only a single comprehensive study of such properties (Van Liedekerke et al., 2010). We have performed a couple of experiments (to be reported elsewhere), in which we have measured the dependence of the length of Epipactis sp. leaves on the applied elongating force. We found that the assumption of the validity of Hooke’s Law (i.e., the linear dependence) is satisfactorily fulfilled; its breakdown happens if the applied forces are close to those which destroy the leaves. The results of our experiments encouraged us to consider the model of a leaf based on a simple version of the theory of elasticity. Our main assumptions in the building of the model are the following: (i) the shape of a leaf is determined by the distribution of its veins, (ii) each vein, together with its surrounding tissue, can be represented by an elastic beam, the shape of which is given by the non-linear theory of strongly deformed beams with circular cross-sections (Landau & Lifshitz, 1993), (iii) the veins are elastically coupled with their nearest neighbours due to the presence of the tissue between the veins, and (iv) only the main (“first-order”) veins are taken into account explicitly whilst the presence of secondary veins may lead to additional concentrated forces acting upon the principal veins. According to our hypothesis, the characteristic undulation near the edges of a leaf is a result of dislocations both in the regions between the veins and in both primary and secondary veins.

In our opinion, the above model has the following justification from the point of view of the leaves micromorphology. In the paper (Jakubska-Busse & Gola, 2014) it has been demonstrated that “due to changes in the number of cell rows (by cell division) and cell sizes/volumes (by cell growth), repeated along the edge of the leaf, sectors with local expansion of the surface affected the entire leaf blade structure” (Jakubska-Busse & Gola, 2014). The resulting “local expansion” can and actually should be physically interpreted as “a defect” in the structure of the beams in the corresponding model. This is because the growth in the number of cells leads to the clear perturbation of the cell sequences—that is, additional, quasi-parallel cell sequences appear rather abruptly. Similar effects appear on the veins themselves, please see additional discussion in ‘Discussion.’

Let dl_m be the infinitesimal element of the arc along the mth beam, and let t_m be the unit vector tangent to that beam. The shape of the beam is given by the function r(l_m) such that ${\displaystyle \frac{d\mathbf{r}}{d{l}_m}={\mathbf{t}}_m.}$

Let F_m denote the force which characterizes the internal stress in the beam. If K_m denotes the external forces (per unit length) which act upon the beam, F has to satisfy (Landau & Lifshitz, 1993): (1)${\displaystyle \frac{d{\mathbf{F}}_m}{d{l}_m}=-{\mathbf{K}}_m.}$

On the other hand, the rate of change of the moment M_m associated with the force F_m along the beam can be written as: (2)${\displaystyle \frac{d{\mathbf{M}}_m}{d{l}_m}={\mathbf{F}}_m\times {\mathbf{t}}_m.}$

For the beam of circular cross section the moment M_m can be written as: (3)${\displaystyle {\mathbf{M}}_m=E_m{I}_m\left({\mathbf{t}}_m\times \frac{d{\mathbf{t}}_m}{d{l}_m}\right),}$ where E_m is the Young modulus of the mth beam, and I_m(z) is its area moment of inertia.

It follows that the tangent vector changes along the beam according to: (4)${\displaystyle E_m{I}_m\left({\mathbf{t}}_m\times \frac{d^2{\mathbf{t}}_m}{d{l}_m^2}\right)={\mathbf{F}}_m\times {\mathbf{t}}_m.}$

Let us obtain the equations for the components of the vector t_m in the spherical coordinates, that is, let ${\displaystyle {\mathbf{t}}_m=\left(cos\left({\theta }_m\right)cos\left({\varphi }_m\right),cos\left({\theta }_m\right)sin\left({\varphi }_m\right),sin\left({\theta }_m\right)\right).}$

Computation of the cross products leads to (please see the Maxima file contained in Supplemental Information 1 for derivation): (5)${\displaystyle E_m{I}_m\left(cos\left({\theta }_m\right)\frac{d^2{\varphi }_m}{d{l}_m^2}-2sin\left({\theta }_m\right)\frac{d{\theta }_m}{d{l}_m}\frac{d{\varphi }_m}{d{l}_m}\right)=F_{x,m}sin\left({\varphi }_m\right)-F_{y,m}cos\left({\varphi }_m\right)}$ and (6)${\displaystyle E_m{I}_m\left(\frac{d^2{\theta }_m}{d{l}_m^2}+\frac{1}{2}sin\left(2{\theta }_m\right){\left(\frac{d{\varphi }_m}{d{l}_m}\right)}^2\right)=F_{x,m}sin\left({\theta }_m\right)cos\left({\varphi }_m\right)+F_{y,m}sin\left({\theta }_m\right)sin\left({\varphi }_m\right)-F_{z,m}cos\left({\theta }_m\right).}$

The above equations have to be supplemented with the the relation: (7)${\displaystyle \frac{d\mathbf{r}}{d{l}_m}={\mathbf{t}}_m,}$ and the expression for the components of the force F_m: (8)${\displaystyle \frac{d{F}_{i,m}}{d{l}_m}=-K_{i,m},}$ where i = x, y, z. We can now define a new dimensionless parameter s such that dl_m = L_ms with L_m being the length of the mth beam. In the following we also use rescaled (dimensionless) forces G_i,m such that: ${\displaystyle F_{i,m}=\frac{E_m{I}_m}{L_m^2}{G}_{i,m}.}$

It is also convenient to work with dimensionless components of the coordinate vector: (9)${\displaystyle {\xi }_m=x_m∕L_m,{\eta }_m=y_m∕L_m,{\zeta }_m=z_m∕L_m.}$

Then the shape of the mth beam together with internal forces are completely specified if we solve the system: (10)${\displaystyle cos\left({\theta }_m\right)\frac{d^2{\varphi }_m}{d{s}^2}-2sin\left({\theta }_m\right)\frac{d{\theta }_m}{ds}\frac{d{\varphi }_m}{ds}=G_{x,m}sin\left({\varphi }_m\right)-G_{y,m}cos\left({\varphi }_m\right),}$ (11)${\displaystyle \frac{d^2{\theta }_m}{d{s}^2}+\frac{1}{2}sin\left(2{\theta }_m\right){\left(\frac{d{\varphi }_m}{ds}\right)}^2=G_{x,m}sin\left({\theta }_m\right)cos\left({\varphi }_m\right)+G_{y,m}sin\left({\theta }_m\right)sin\left({\varphi }_m\right)-G_{z,m}cos\left({\theta }_m\right).}$

From the definition of the tangent vector t_m we also have: (12)${\displaystyle \frac{d{\xi }_m}{ds}=cos\left({\theta }_m\right)cos\left({\varphi }_m\right),}$ (13)${\displaystyle \frac{d{\eta }_m}{ds}=cos\left({\theta }_m\right)sin\left({\varphi }_m\right),}$ (14)${\displaystyle \frac{d{\zeta }_m}{ds}=sin\left({\theta }_m\right).}$

We have assumed the following simple form of the force ${\bar{\mathbf{K}}}_m=L_m^3∕\left(E_m{I}_m\right){\mathbf{K}}_m$: (15)${\displaystyle {\bar{K}}_{i,m}=-{\kappa }_{i,m}\left(2{\rho }_{i,m}-{\rho }_{i,m+1}-{\rho }_{i,m-1}\right)-q_{i,m}{\delta }_{i,z},}$ where i = x, y, z and q_i,m = q_mδ_iz is the scaled, dimensionless weight of the mth beam per unit length. Therefore, the equations for the dimensionless forces G_i,m take the form: (16)${\displaystyle \frac{d{G}_{i,m}}{ds}={\kappa }_{i,m}\left(2{\rho }_{i,m}-{\rho }_{i,m+1}-{\rho }_{m-1}\right)-q_m{\delta }_{i,z}.}$

In the following we have assumed that there is no external force with non-vanishing x coordinate acting on any beam except, perhaps, of that acting at the tip of the leaf. Thus, the above system of differential equations simplifies a bit since G_x,m is independent of s. The resulting system of differential equations has the order 9M and its solution is available only numerically. The boundary conditions has been chosen as follows. We have used η_m(0) = ζ_m(0) = 0, ξ_m(1) = 1, η_m(1) = ζ_m(1) = 0, G_y,m(1) = G_z,m(1) = 0. We have also prescribed ϕ_m(0) = ϕ_m0, θ_m(0) = θ_m0 and experimented with various ϕ_m0 and θ_m0. The above boundary conditions correspond to the left (s = 0) end of each beam lies on the ξ-axis. About its right end (s = 1) one can say that the sum of all forces acting on the beams there vanishes. As the tip can be considered point-like, the beams meet at one point. We can assume (performing rotation of the coordinate system if necessary) that the tip is located at (1, 0, 0), hence the boundary conditions for ξ_m(1), η_m(1), ζ_m(1). The above boundary condition effectively imply that we have to do with the buckling of beams in our model due to the concentrated forces very close to the tip of the leaf. That buckling clearly influences, to some extent, the shape of the blade as well as the wrinkling. Ideally, we should have prescribed a (possibly much more complicated) correct model of the inter-beam forces. The convergence of the beams near the tip would then be a consequence of such a good model. Unfortunately, we have encountered serious difficulties in finding such a convincing description. So far, we have decided to introduce the tip simply by the above boundary conditions. It seems to us, however, that the presence of strong concentrated forces near the end of the blade is not unlikely, and that the resulting buckling does influence the shape of the leaf. Let us notice that the ξ coordinates of the left ends of the beams are unspecified and they actually result from the simulation. They are always close to, but never exactly equal to 0. This results in appearance of a small “petiole.”

In addition to forces specified above, we have assumed the presence of concentrated, Dirac-delta like forces in order to model dislocations (this is consistent with the elementary theory of dislocations as described in Landau & Lifshitz (1993)). That is, we have added the forces of the form ${\displaystyle f_{z,m}\left(s\right)=\sum _k{g}_{m,k}\delta \left(s-s_k\right)}$ to the right-hand sides of Eq. (16). Only the z-component of (16) has been modified this way.

Apart from gravity, the total external force which acts on the beams must naturally be equal to zero.