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

Static properties of leaves with parallel venation, with particular emphasis on the genus EpipactisZinn, 1757 (Orchidaceae, Neottieae) have been modelled with coupled quasi-parallel elastic “beams.” The non-linear theory of strongly bended beams have been employed. The resulting boundary-value problem has been solved numerically with the help of the finite-difference method. Possible dislocations resulting in additional Dirac-delta like forces have been take into account. Morphological similarity of the model and real leaves has been obtained. In particular, the concentrated forces have been shown to cause undulation in the leaves.


INTRODUCTION
The modelling of plant organs remains an open problem due to the complexity of plant architecture, regardless of the particular organ or any particular taxonomical group of plants being considered. In this study, we are concerned with the leaves of Orchidaceae and specifically, the leaves of plants belonging to the genus Epipactis.
Orchidaceae are a varied group that includes both terrestrial plants and epiphytes with monopodial and sympodial growth, which differ from each other by leaf construction in both shape and size. Most orchid leaves are typical of the monocots with many parallel veins and the cross connections between the longitudinal veins being inconspicuous (Dressler, 1981;Jakubska-Busse & Gola, 2014). An orchid leaf usually consists of only a blade, though some species produce also the petiole which may be narrow or which may widen and embrace a stem, forming a leaf sheath. The shape of the blade tends to be modified alongside its size as a result of growth processes (Jakubska-Busse & Gola, 2014). Orchid taxa may produce leaves of different shapes, e.g., plicate and cordate, plicate with a sheathing base, convolute, non-articulate, conduplicate and deeply lobed or, rarely, twisted (Dressler, 1981).
An interesting and relatively well-studied issue is leaf undulation, which is often observed in monocots including orchids. It has been found that the waves in the leaf blades in monocots usually appear perpendicular to the leaf length, which demonstrates that as the leaf surface grows, it changes correspondingly lengthwise to the leaf blade. It is worth noting that although undulation normally occurs perpendicularly to the leaf length, at the beginning stage it occurs alongside its length. Displacements in the wrinkled leaf also occur across the leaf blade. Whichever way they appear, they demonstrate that the pace of their growth is irregular since wrinkles in a leaf can be of different length. Hejnowicz (1991) found that spatial and temporal fluctuations in the pH of the epidermal cell walls aided the undulation.
Epipactis as the main object of our study is a mainly Eurasian genus with a southcentral distribution (Rasmussen, 1995;Delforge, 2006). These are all rhizomatous (clonal), summer-green plants with habitats ranging from bogs to dry forest (Delforge, 2006). There are some controversies around the taxonomy as one of the Epipactis muelleri Godf. species is identified based on the undulation of its leaf blade, alongside gynostemium features. This property is quoted in nearly every artificial key for determining this group of orchids, e.g., Sundermann (1975), Procházka & Velísek (1983), Potůček & Čačko (1996 and Delforge (2006). However, it has been shown by Jakubska-Busse & Gola (2014) that the leaf undulation in Helleborines does not have any diagnostic value as a non-programmed intrinsic feature and should not be applied to taxa identification. In this paper, we attempt to show which characteristics in the mathematical model of leaf are responsible for the undulation.

MATERIAL AND METHODS
In order to build the model, we have preliminarily checked the fulfillment of Hooke's law on a sample of 72 fresh Epipactis sp. (E. helleborine, E. muelleri, E. albensis and E. palustris) leaves. Those measurements have been used here to justify the approach based on theory of elasticity as modeling tool. (We plan to perform much more extensive measurements in some future.) A code written in the Python language was used to solve the model described below and this employed the finite-difference method to solve boundary-value problem for a system of ordinary differential equations. We have also used published data on the morphology and anatomy of the Epipactis sp. leaves (Jakubska-Busse & Gola, 2010;Jakubska-Busse & Gola, 2014;Jakubska-Busse et al., 2012). The presented studies were done with the consent of the Regional Director for Environmental Protection, No. WPN.6400.27.2015.IW.1.

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 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): 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: For the beam of circular cross section the moment M m can be written as: 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: Let us obtain the equations for the components of the vector t m in the spherical coordinates, that is, let t m = (cos(θ m )cos(φ m ),cos(θ m )sin(φ m ),sin(θ m )).
Computation of the cross products leads to (please see the Maxima file contained in Supplemental Information 1 for derivation): and The above equations have to be supplemented with the the relation: and the expression for the components of the force F m : where i = x,y,z. We can now define a new dimensionless parameter s such that dl m = L m s with L m being the length of the mth beam. In the following we also use rescaled (dimensionless) forces G i,m such that: It is also convenient to work with dimensionless components of the coordinate vector: Then the shape of the mth beam together with internal forces are completely specified if we solve the system: From the definition of the tangent vector t m we also have: We have assumed the following simple form of the forceK m = L 3 m /(E m I m )K m : 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: 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  (1993)). That is, we have added the forces of the form 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.

NUMERICAL RESULTS
The above boundary-value problem has been solved using its representation in terms of four second-order equations and finite differences. The discretization has been performed with N = 101 points and s i − s i−1 =: h = 1/(N − 1). It has been assumed as a kind of ''zeroth-order approximation'' that all q m are the same, q m = q = −1.0.
The coefficients g (m,k) have been assumed on m as follows: so that the (scaled) internal concentrated forces are larger near the edges and smaller in the central part of the leaf. Since these forces are scaled, that difference can be attributed either to forces themselves, or to smaller Young modulus and/or area inertial moments of the lateral beams. We have worked with M = 32 beams and used the finite differences to discretize the differential operators. Since the resulting matrices in our boundary-value problem have been quite large, we have not used any solver, but employed the method of false transients instead. The results of these simulations are illustrated in Figs. 1-4 for four sets   Figure 1 displays the results for a leaf under the absence of any defects (hence, no point-like forces appear) and small coupling between the beams which form the leaf blade.
In Fig. 2 we have shown a shape which has resulted from much stronger coupling between the beams, but still under the absence of any point-like forces. The difference between the shapes in Figs. 1 and 2 is evident, but there is no undulation.
We have performed preliminary measurements of the amplitudes of wrinkles for several leaves. It has been found that the amplitudes of the leaf wrinkles have normally been of the order of 1-2% (it's depended of course on whether we measured in the apical, central or basal part of the leaf) but not exceeded 10% of its length The parameters of the model are chosen in such a way as to fit that maximum value of the relative (i.e., related to the length of the leaf) amplitudes. In Fig. 3, however, the number of defects as well as the strengths of point-like forces are too small to cause any considerable wrinkles. Only in Fig. 4 are they strong enough to produce a clear picture of undulation. As we observed before, large values  (i.e., non-vanishing g 0 ). The following values of parameters have been used: q 0 = −1.0, θ m (0) = 0, m = 1,2,...,M , κ y = 10.0, κ z = 10.0, s 2i−1 = i/10,s 2i = (i + 1/2)/10,i = 1,2,...,6; g 0 (1) = g 0 (3) = −g 0 (2) = −g 0 (4) = 500,g 0 (5) = g 0 (7) = g 0 (9) = g 0 (11) = −g 0 (6) = −g 0 (8) = −g 0 (10) = −g 0 (12) = 400. The parameters correspond to strong defects or a small Young's modulus close to the edges, and it results in stronger undulations in those regions. of the parameters g 0 (i) can be attributed to both strong defects and to small values of the Young modulus of the beams. In particular, smaller values of the Young modulus close to edges results in the stronger undulations in those regions. Please see Fig. 5 for comparison.

DISCUSSION
In this paper the model of a plant leaf with quasi-parallel venation (as exemplified by the Epipactis sp. leaf) which consists of coupled elastic beams has been developed. Only the nearest-neighbour interactions between the beams have been employed. It was found that the model reproduces quite well the overall shape of the leaf. In particular, it reflects to some extent the presence of undulation at the leaf edges. The appearance of undulation is attributed to the following factors: (a) a small Young's modulus and/or moment of inertia of veins which are close to the edges; (b) dislocations in the structure of the principal veins themselves and/or connecting tissue and secondary veins which result in the large but well-localized additional forces which act on the veins; those dislocation appear during the growth of the leaf.
The presence of the defects in the sequences of cells in the space between veins has been discussed in (Jakubska-Busse & Gola, 2014). As for the presence of clear change of such sequences close to the veins themselves, let us just invoke Fig. 6 where it seems to be transparent.
Some of the shortcomings of the analysis in this paper are the following. Firstly, even though the theory of dislocations in elastic bodies predicts the presence of Dirac-delta-like forces, it is entirely unclear whether they could act in precisely the way we have described above. Detailed microscopic observations together with thorough and complex micromechanical models could bring the answer. Quite obviously, we have applied the Dirac-delta-like forces in an ad hoc manner to obtain numerically results close to those known in real leaves.
Secondly, we have assumed a lot (and possibly too much) about our beams to simplify the model. Needless to say, the cell sequences in the leaves do not have circular cross sections. It may very well happen that this fact, together with strong forces acting near the apex of a leaf, is a source of instability which is and additional factor in formation of the undulation. Nonetheless, we believe that our arguments based on the findings of Jakubska-Busse & Gola (2014) and Fig. 6 provide first necessary steps in the proper direction.
In some remarkable recent works (Marder, 2003;Audoly & Boudaoud, 2004;Sharon, Roman & Swinney, 2007), the variety of shapes of leaves in general, and their wrinkling in particular, has been interpreted in a very elegant way, namely, in terms of Riemannian geometry of surfaces. What is more, those authors have described the effect of applying auxin at leaf edges to produce ripples and emphasize more explicitly growth in the leaf. It is conjectured that in the developmental processes of naturally growing tissues, the process of growth provides a mechanism for the spontaneous formation of non-Euclidean metrics and consequently leads to complicated morphogenesis of thin films (including surfaces of leaves) exhibiting waves, ruffles and non-zero residual stress (Lewicka, 2011). What we have attempted here is to provide a more physical, granular model which is from the very beginning adapted to the strong anisotropy of Epipactis leaves and concerned with local forces and torques. We believe that, to some extent, our model is complementary to that developed in the papers quoted above, being, perhaps, somewhat closer to the microscopic reality of the Epipactis leaves. An interesting question which arises in this connection is whether we could appropriately interpret within our model the results of the experimentation with auxins. Let us first stress that the leaves used in the those experiments (daffodils (Narcissus)) are quite different from those of Epipactis. Nevertheless, our model seems to be sufficiently abstract to justify its application to other leaves with both similarly shaped leaves and similar arrangements of main veins. The outcomes of experiments with auxins could be explained within our model as resulting from the (local) change of physical characteristics of the beams, especially their area moments of inertia. We also feel that the departure from the assumption of the circularity of cross sections of the beams might be necessary.
Further research is necessary in at least two directions. Firstly, the experimental data regarding the mechanical properties of the leaves in general, and the leaves of Orchidaceae in particular, is presently insufficient. Secondly, it is obviously necessary to combine a purely mechanical model with hydrodynamics and the diffusive properties of the dynamics of fluids inside the leaf. We plan to continue our studies in both of these directions.
• Maciej Janowicz conceived and designed the experiments, analyzed the data, wrote the paper, prepared figures and/or tables, reviewed drafts of the paper.
• Luiza Ochnio and Beata Jackowska-Zduniak conceived and designed the experiments, performed the experiments, analyzed the data, contributed reagents/materials/analysis tools, wrote the paper, reviewed drafts of the paper.

Ethics
The following information was supplied relating to ethical approvals (i.e., approving body and any reference numbers): The presented studies were done with the consent of the Regional Director for Environmental Protection, No. WPN.6400.27.2015.IW.1.

Data Availability
The following information was supplied regarding data availability: The raw data has been supplied as Supplemental Information 1.

Supplemental Information
Supplemental information for this article can be found online at http://dx.doi.org/10.7717/ peerj.2165#supplemental-information.