A Multiscale Mechanical Model for Plant Tissue Stiffness

Plant petioles and stems are hierarchical cellular structures, displaying structuralfeatures defined at multiple length scales. The current work focuses on the multi-scalemodelling of plant tissue, considering two orders of structural hierarchy, cell wall and tissue.The stiffness of plant tissue is largely governed by the geometry of the tissue cells, thecomposition of the cell wall and the structural properties of its constituents. The cell wallis analogous to a fiber reinforced composite, where the cellulose microfibril (CMF) is theload bearing component. For multilayered cell wall, the microfibril angle (MFA) in themiddle layer of the secondary cell wall (S2 layer) largely affects the longitudinal stiffnessfor values up to 40o. The MFA in turn influences the overall wall stiffness. In this work,the effective stiffness of a model system based on collenchyma cell wall of a dicotyledonousplant, the Rheum rhabarbarum, is computed considering generic MFA and volume fractions.At the cellular level, a 2-D Finite Edge Centroidal Voronoi tessellation (FECVT) has beendeveloped and implemented to generate the non-periodic microstructure of the plant tissue.The effective elastic properties of the cellular tissue are obtained through finite elementanalysis (FEA) of the Voronoi model coupled with the cell wall properties. The stiffness ofthe hierarchically modeled tissue is critically important in determining the overall structuralproperties of plant petioles and stems.


Introduction
Plants are one of the major kingdoms in biology.From a structural point of view, a plant exhibits remarkable mechanical properties.Among their organs, the petiole that attaches the leaf to the stem is one of the significant load bearing structures of a plant.As a cantilever beam, the petiole supports the leaf against gravity and allows it to be exposed to the sun.The petiole thus provides mechanical support against the weight of the leaf and against environmental factors, such as rain and wind, resisting both bending and twisting load [1].The petiole of Rheum rhabarbarum (rhubarb) plant, shown in Figure 1, is an example of a cantilever beam that must resist combined loads including bending and twisting.The petiole's flexural and torsional stiffness are largely influenced by its overall geometric properties and the stiffness of its constituents, including tissues and cell walls.The petiole can be considered as a hierarchical structure, having structural features defined at multiple length scales.A petiole is generally composed of an assembly of cellular tissues whose mechanical characteristics collectively depend on the geometry of their constituent cells, cell wall composition, the structural properties of the wall constituents, and the microstructures of the tissues.It is useful to characterize the structural hierarchies in plants in terms of levels (n).As indicated in Figure 2, n = 1 corresponds to the cell wall, n = 2 to tissue, and n = 3 to the petiole.Although the hierarchical organization of plants spans from the molecular level to the whole organism, in this work the stiffness of tissues is specified by the hierarchies at the cellular and sub-cellular level (Figure 2).Thus, the mechanical attributes of a tissue depend on the structural properties of the cell wall constituents and their structural organizations.
Cell walls are mainly composed of complex networks of polysaccharides, namely cellulose (C), hemicelluloses (HC) and pectin, along with comparatively minor quantities of structural protein and/or lignin [2][3][4][5].Each Cellulose Microfibril (CMF) is a semi-crystalline aggregate of 36 β-1,4-glucan chains.CMFs are several µm long but only about 3-4 nm thick, and lay 15-30 µm apart in the wall [6].The construction of the cell wall is considered to be a fiber-reinforced composite, where the cellulose fibrils act as the main load-bearing elements [7][8][9].The CMFs are embedded in a compliant matrix of hemicelluloses and lignin, where the stiffness of the matrix is approximately two orders of magnitude lower than the CMF [4,10].Experiments on plant cells show that the properties of CMF strongly dictate the elasticity of the cell wall, while lignin and hemicelluloses have a marginal effect [11].However, the stiffness of cell wall can be approximated considering the simplest composite model, where the cell wall is a two phase system under uniaxial loading condition [12,13].Based on the elastic properties and the volume fractions of the CMF and lignin/hemicellulose matrix, the stiffness of the cell wall can be determined for two limiting architectures-parallel (Voigt model) and perpendicular fiber orientations (Reuss model)-with respect to loading direction.Since there is a significant difference in stiffness between the polysaccharides at the cell wall level and other levels of hierarchy, the orientation of the microfibrils in the cell wall influences the stiffness properties of individual cells and tissues as well.The Microfibril Angle (MFA) has been experimentally shown to be a key factor in determining the mechanical properties of plant cell walls [3][4][5]9].Both experimental and numerical studies show that the structural properties of plants are influenced by the properties of cell wall constituents, MFA and their structural organization [14][15][16].Hence, the Voigt and Reuss models are insufficient to approximate the stiffness of the cell wall.To determine the longitudinal and transverse modulus, the orthotropic nature of the cell wall should be considered [17,18].The stiffness of the cell wall should thus be modeled assuming orthotropic material properties, where the constitutive relation yields the stiffness of the cell wall along the orthogonal directions.It has been demonstrated that the shape, size, and spatial distribution of cells govern the physical, biological, and structural properties of cellular materials [20][21][22].Hence, the ability to realistically model the cellular microstructure of a plant tissue is crucial to understanding its mechanical behaviour [23,24].Many researchers have modeled natural cellular solids using repeating unit cells to construct a regular microstructure in the form of a circular, square or hexagonal array of cells [25].Closed-form relations of the structure-property can be derived with simplified geometric models based on repeated unit cells [26,27].Nevertheless, the modeling of plant tissue is a challenging task, since natural cellular solids often exhibit non-periodic arrangement of cells.Since the microstructure is heterogeneous in shape and size, the Voronoi tessellation can be used to generate an accurate representation of a non-periodic microstructure [27][28][29].Voronoi tessellation is used extensively to model grain geometry for the property characterization of polycrystalline aggregates [30] and inter-granular cracks [31].Mattea et al. [32] and Roudot et al. [33] pioneered the use of Voronoi tessellation to model the microstructure of fruit tissues.Both groups aimed only to capture the randomness of the fruit tissue microstructure without necessarily producing a model that accurately represented the real tissue.Recently, Ntenga et al. [34,35] tried to analyse the structure, morphology and mechanical properties of the Rhectophyllum camerunense (RC) plant fiber using a conventional Voronoi diagram.Due to inherent drawbacks of the Voronoi (CVT) model, semi-infinite edges were present at the boundary of the fiber, making the model unsuitable for Finite Element Analysis (FEA).To generate a geometric model, having finite edges at the boundary, a Finite-edge Centroidal Voronoi Tessellation (FECVT) method has been developed.The FECVT method has been applied to model the tissue microstructures of emphArabidopsis thaliana and Philodendron melinonii, two plants with distinct tissue architecture [36].In this work, the FECVT method is applied to generate the virtual model of the rhubarb tissue.Coupled with the cell wall properties, the FEA of the virtual model provides the effective/homogenized tissue stiffness.
Most studies on the homogenization of cellular solids are based on regular models with a periodic microstructure.Real solid foams, however, exhibit amorphous arrangements of pores with different sizes and shapes rather than perfectly periodic structures.The homogenized/apparent elastic property of periodic honeycomb varies from 10% to 15% compared with non-periodic honeycomb [27].On the other hand, to analyze the irregular microstructures represented by the Voronoi tessellation, the homogenization process requires a large scale representative volume element (RVE).Analyses of such models had been provided, among others, by Silva et al. [27], Fazekas et al. [37], Roberts and Garboczi [38] for both two and three-dimensional models.To obtain the homogenized/effective stiffness properties, the 2-D virtual tissue (RVE) undergoes FEA for a given boundary condition [39].The computational homogenization technique is used with non-periodic microstructure for global homogenization.This technique can effectively consider the microstructural irregularity present in the plant tissues.The effective properties can be used in the next hierarchical level to capture macro-scale properties such as bending and torsional stiffness of the whole petiole.
The objective of this work is to determine the theoretical bounds of stiffness of the collenchyma tissue modeled with two orders of structural hierarchy, namely cell wall and tissue micro-architecture.At the first order of the hierarchy, the constitutive properties of the cell wall are modeled via the classical micromechanical models of Voigt and Reuss, and the microfibril angle.The geometric model of the collenchyma tissue is generated by the FECVT algorithm; it undergoes finite element analysis that accounts for the cell wall properties to determine the stiffness of the tissue at the second order of the hierarchy.

Methodology
Since most plant tissues display hierarchical organization, we adopt in this work a hierarchical modeling approach to obtain their stiffness properties.In the first order of the hierarchy, the cell wall stiffness is modeled considering as variables the MFA and the volume fractions of the wall constituents, and their elastic properties.The obtained wall stiffness is then used in the next hierarchical level, where the effective stiffness of tissue, represented by an FECVT model, is determined through FEA.

Stiffness of Cell Wall via Micromechanics Theory
The cell wall, which belong to the bottommost hierarchy (n = 1 in Figure 2), is the building block at this level of hierarchy.The upper and lower limits of stiffness of the cell wall can be determined using micromechanical models-Voigt and Reuss models.Three pieces of information are required in the model: (i) the elastic moduli of cell wall constituents; (ii) the volume fractions of the various constituents; and (iii) the geometric arrangements of the constituents relative to each other.The properties of the cell wall constituents are available in literature.For example for cellulose, a wide range of theoretical as well as experimental estimation of its elastic modulus is reported.Earlier theoretical prediction of Young's modulus provides relatively higher values, 246 and 319 GPa [40,41], whereas more recent analyses show 120-140 GPa [42][43][44] and 167 GPa [45] for crystalline cellulose along the fiber direction.Data for the elastic modulus of hemicelluloses are scarce and fall within a range of 5-8 GPa [46].The elastic properties of lignin are difficult to assess and vary from 4 to 7 GPa [47,48].Considering an orthotropic arrangement in the cell wall, a set of rational moduli [10] used in this work are shown in Table 1.The volume fractions and the geometric arrangements of the constituents are discussed in Section 3.1, where their effects on the engineering properties of cell wall are discussed in detail.
Table 1.Elastic properties of cell wall constituents used in the present work [10].

CMF
Lignin Hemicellulose The values shown in Table 1 are used in the Voigt and Reuss models to determine the wall stiffness along the fiber direction, E 1 , and across the fiber direction, E 2 (Figure 3).The Voigt upper bound is and the Reuss lower bound is where E 1 and E 2 refer to the Young's moduli of the composite, E f and E m refer to the Young's modulus of fiber and matrix, respectively; v f and v m refer to the volume fraction of fiber and matrix, respectively, with Using the same concept, the major Poisson's ratio, ν 12 , along the fiber direction, is and the minor Poisson's ratio, ν 21 , across the fiber direction, can be derived with the reciprocal relationship The in-plane shear modulus can be expressed as However, an actual cell wall may not behave according to the predictions based on either the pure Reuss or the pure Voigt model.Rather, a hybrid model appears to be the best modeling approach although the Voigt elements appear to be dominant.

Constitutive Properties of Cell Wall using Macromechanics
The constitutive relation for the CMF (fiber) reinforced cell wall depicts the stiffness in the global coordinate system.For an orthotropic cell wall, the compliance matrix is expressed for plane stress condition as where The constitutive relation for an orthotropic material is therefore, where Equation ( 7) stands for the on-axis configuration, where the global coordinate system, x-y axis, coincides with the local coordinate system, 1-2 axis, respectively.When the local and global coordinates do not coincide, the arrangement is the off-axis configuration, shown in Figure 4. Therefore, the on-axis stress-strain relationship is not adequate for the analysis of an off-axis configuration.In the off-axis configuration, the stiffness depends on the fiber orientation angle, which is termed as the MFA for the cell wall configuration.Thus, the MFA needs to be addressed in the constitutive relation.If the x-y coordinate system is aligned with the direction of loading and the 1-2 coordinate system is aligned with the fibers, MFA (θ) is the angle between the two (Figure 4).
The local (1-2) and global (x-y) coordinate systems are related through a transformation equation.Hence, the transformation of strains between the loading (local) direction and the fiber (global) direction is denoted by where m = cosθ, m = sinθ and [T ] is the strain transformation matrix.Similarly, the stress transformation matrix is Reorganizing and combining Equations ( 7)-( 9), the Cartesian strain-stress relationship in the global coordinate system is given by where [ S] is a new matrix known as the transformed reduced compliance matrix.The off-axis constitutive relationship can be expressed in the following form: In summary, the constitutive (engineering) properties of cell wall are

Geometric Modeling
A Voronoi microstructure is constructed based on a set of randomly generated points called Voronoi sites.The cell boundaries are drawn such that any point within the enclosed polygon is closer to its Voronoi site than to the Voronoi sites of the surrounding polygons.The Voronoi tessellation thus divides a space into as many regions as there are Voronoi sites.In the current method, the centroids of the cells are used to construct the Voronoi diagram.To construct a centroidal Voronoi tessellation (CVT), first, a color or greyscale micrograph of a plant tissue is subjected to image segmentation.Thresholding, a method for image segmentation, converts the color or greyscale image into a binary (black and white) image.Since plant tissues usually exhibit graded cellularity as well as complex heterogeneity, an edge detection algorithm is used in conjunction with thresholding to obtain the cellular distribution accurately.The Canny edge detection algorithm [49] is used here to detect the shapes of the cells more precisely.
Once the shapes of the cells are detected, the first order moments of the cells are computed using X and Y coordinates of the pixels.The algebraic form of the moment equation of a digital image is where (X i , Y j ) is the coordinate of the (i, j) th pixel, f (i, j) have value 1 if the (i, j) th pixel is in the shape and 0 otherwise.Considering the region of interest, which is completely enclosed in a rectangular region G of size n 1 by n 2 pixels, i varies from 1 to n 1 , and j varies from 1 to n 2 in the function f (i, j).
For a 2-D region, p + q denotes the order of moment, where p and q are integers.Hence, the coordinates of a cell are X = m 10 m 00 and Ȳ = m 01 m 00 (18) where the zeroth moment, physically, is equal to the area of the region.
After determining the centroids of the cells, the Voronoi tessellation is constructed based on the Quick hull algorithm [50].The outcome is a conventional CVT with semi-infinite edges at the boundary.The semi-infinite edges complicate finite element analysis because the boundary conditions, applied at an infinite distance, are not realistic.This problem is especially difficult to correct in models with irregular shape contour.
To remove the infinite edges from the boundary, the centroids of the outermost cells are determined.For each centroid, the distances between the centroid and the surrounding Voronoi sites (centroids of the surrounding polygons) are calculated, and the minimum distance is determined.An imaginary point is created such that the distance between itself and the centroid is half of the minimum distance.The imaginary point is thus created for each of the selected centroid.Based on the set of imaginary points and the convex hull algorithm [50], a boundary is imposed, after which a Boolean subtraction is realized.With this Boolean operation, the semi-infinite edges are truncated, and the vertices of the truncated edges are reconnected to form the final boundary.The FECVT technique is, therefore, capable of capturing the microstructure of an image with an arbitrarily shaped boundary contour.This novel computational algorithm, FECVT method [36], is applied here to model the rhubarb tissue microstructure.

Finite Element (FE) Modeling
The compressive stress-strain regime of the Voronoi microstructure is determined using finite element analysis (ANSYS Inc., Canonsburg, PA, USA).Each cell wall of the Voronoi microstructure is modeled with a BEAM23 element, capable of modeling both elastic and plastic behaviour.The shear deformation, which is important for stubby beam, is also captured by considering the shear deflection coefficient of the beam element.Beams have a rectangular cross-section of uniform thickness, t.The relative density ρ * /ρ s of the microstructure is given by where N is the total number of beams, l i is the length of the beam i; L X and L Y are the dimensions of the Voronoi model along X and Y directions, respectively.The FEA is conducted for different relative density, adjusted by the value of t.In the finite element analysis, a Young's modulus E s = 1 is assigned to each beam to obtain the normalized tissue stiffness.

Boundary Conditions
Three types of boundary conditions (BC)-(1) mixed boundary condition, (2) prescribed displacement boundary condition, and (3) periodic boundary condition-can be generally imposed.Since the microstructure is not periodic, the periodic boundary condition is not appropriate here.The mixed boundary condition enforces the normal displacement, eliminating the tangential force and the bending moments at the nodes on the boundaries.Since the mixed boundary condition tends to underestimate the Young's modulus [51], the displacement boundary condition has been selected, which has been proved to be appropriate for the analysis of non-periodic microstructures for uniaxial and biaxial loading cases [27,29,37,52].
To determine the effective Young's modulus, uniaxial compression load with prescribed displacement is simulated.A uniaxial compressive strain along the X axis is imposed on the nodes of the right edge while the nodes at the left edge are constrained from translating in this direction (Figure 5).The nodes of the bottom edge are also constrained from translating in the Y direction to prevent the rigid body motion.Similarly, uniaxial compression along Y axis has also been performed.In both cases, nodes are constrained from rotating in the X − Y plane.To determine the effective shear modulus, a biaxial loading test is simulated with a positive displacement in the X direction and a negative displacement in Y direction.The results are computed for alternative values of densities for prescribed boundary conditions.

Determination of Apparent Tissue Stiffness
For each model with a given relative density, the effective Young's modulus, E * , and effective shear modulus, G * XY , are determined.The macroscopic stress, σ * , is calculated from the global reaction forces of the structure in the loading direction.The sum of the nodal reaction forces is divided by the edge length to determine the normal stress in the loading direction.The strain, * , is determined based on the technique of gage lines introduced by Silva [27] to eliminate the end effects.The shear strain is computed as the change in the angle between the gage lines oriented at 45 • , with respect to the coordinate axes.

Results and Discussion
The multiscale hierarchical approach can be applied effectively to simulate the stiffness of various plant tissues.In the present work, the method is implemented to the collenchyma tissue of a dicotyledonous plant organ, the rhubarb petiole.

Stiffness of Cell Wall
The properties of a cell wall depend on the wall constituents, mainly cellulose microfibril, lignin, and hemicellulose.The CMF of cellulose acts as the load bearing element that controls the wall properties.Along with the elastic properties of the constituents (shown in Table 1), the volume fractions of the constituents play an important role as well.The volume fractions of the constituents are not homogenous throughout a multilayered plant cell wall.The layered construction is usually evident in woody plants but is also observed in many other plants.Although the number of layers depends on several factors, such as the types of species, cells, the cell wall is widely accepted to be considered as a three-layered structure.In fact, the woody plants typically exhibit three secondary layers along with a primary layer and middle lamellae.The primary layer only plays a significant role in the early stage of a growing plant.When a plant cell reaches its maturity, the secondary wall becomes the determinant governing its mechanical performance.The volume fractions of CMF, lignin, and hemicelluloses shown in Table 2 represent values of a typical set of multilayered plant cell walls.However, the composition of the constituents depends on other factors as well, such as types of species, and cells, location of the cells and the age of a plant.The data shown in Tables 1 and 2 are used to compute the effective engineering properties, such as stiffness and Poisson's ratio, of cell wall.The variations of stiffness with respect to MFA along different layers are shown in Figure 6. Figure 6a depicts the effect of MFA on longitudinal the Young's modulus for each layer.Up to 40 • , the effect of MFA on the longitudinal wall stiffness is substantial, and the stiffness is reduced to a considerable amount.In contrast, the transverse stiffness is nearly invariant with the MFA up to 50 • .Beyond this limit, the transverse modulus varies with MFA and increases with the increasing MFA.The shear modulus is maximum for CMF equal to 45 • in the loading direction.However, the variation of shear modulus up to 20 • is less significant, and a gradual increase is observed until the MFA reaches 45 • and starts decreasing afterwards.The Poisson's ratios shown in Figure 7 vary with the MFA.From the comparison of the minor Poisson's ratio, we observe that the major Poisson's ratio increases rapidly up to MFA of 20 • .The major Poisson's ratio plays an important role in influencing the anisotropic cell wall properties.
It is evident that the MFA plays a significant role in the wall stiffness, an observation that has been experimentally noticed in numerous plants' cell wall.The MFA usually varies between 0 • and 40 • , in the S 2 layer [53][54][55] and plays a leading role in determining the overall wall stiffness.Moreover, since the thickness of S 2 layer is generally much higher than that of S 1 and S 3 (also primary layer) layers, the measurement of MFA for the whole cell wall or the average MFA across the cell wall involves the approximation of MFA in the S 2 layer.Thus, the S 2 layer is often the main determinant of the stiffness of plant [56].This concept is used here to determine the effective engineering properties of a single layered dicotyledonous cell.Although a more accurate and rigorous analysis would be conceivable considering microfibril angles at different layers, the composition of S 2 layer with respective MFA provides a good approximation for the overall stiffness of cell wall.Moreover, the measurement of microfibril angles in different layers is not always experimentally possible.Since the focus of this paper is to develop a multiscale model that accounts for the elastic part of the deformation, the nonlinear mechanical response has been ignored.Yet, the non-linear response of plant tissue is important because the approach presented here can be extended to model the nonlinear deformation of the cell wall as well.In stems and petioles, collenchyma tissue is observed in peripheral locations, beneath the epidermal layer (Figure 8).As per experimental measurements, the collenchyma tissue typically comprises 20% cellulose and the rest 80% of HC and lignin [57,58].Although there is no clear distinction, the collenchyma cell walls can be considered to be secondary cell wall [58,59], where the orientation of CMF is mostly found to be inclined to the longitudinal axis [60][61][62][63][64][65].However, the MFA also varies with the technique used for the measurement.For dicotyledonous plants, the MFA usually lies between 6 • and 15 • [66].Based on the fiber volume fraction and MFA, the stiffness of collenchyma tissue is derived.The effective stiffness properties of collenchyma cell wall with 15 • MFA are summarized in Table 3.The effective stiffness of the collenchyma tissue will be determined based on the effective wall stiffness.

FECVT Modeling of the Tissue Microstructure
The geometry of the collenchyma tissue of the rhubarb petiole is generated using the FECVT algorithm.Figure 9 shows the collenchyma tissue and its corresponding FECVT model.The latter is obtained with a method recently introduced to represent complex geometric heterogeneity in plant tissue [36].Since experimental data are used to formulate the mathematical model, the accuracy of the model appears to depend on the quality of the micrograph.The FECVT method, thus, can capture the details of a cellular distribution if the micrograph of the tissue microstructure is vivid and clear.If a conventional CVT were to be used, the reproduction of the polygons at the boundary would differ for their shape and size from those obtained with a FECVT.FECVT models are generated for several arbitrarily chosen sections of collenchyma tissue.Each section is a representative of the respective tissue.Nonetheless, a one-way ANOVA is used to test the differences of FECVT models generated for the tissue and shown in Table 4.The first column in Table 4 corresponds to the between-groups and within-groups estimate of variance, which are shown in the fourth column, showing the mean square.The second column gives the sum of squares for each of the estimates of the variance.The third column gives the degrees of freedom (df) for each estimate of the variance.The number of FECVT models and the sum of the Voronoi cells in the models are used to calculate the degrees of freedom between-groups and within-groups, respectively.The fifth column, F, corresponds to the ratio of the variance of between-groups and within-groups, shown in the preceding column.The last column gives the probability of significance of the F ratio, i.e. p value.Since the p value is greater than the significance level (α = 0.05), the null hypothesis is not rejected.Therefore, the FECVT models statistically do not differ significantly across the models as determined by one-way ANOVA (F (2, 611) = 0.172, p = 0.837 > 0.05).The output of the ANOVA test implies that the FECVT models of the stochastically chosen collenchyma tissue sections are statistically invariant; hence each model can be considered to be a representative volume element for the FEA.Table 3. Stiffness properties of collenchyma cell wall determined via Equations ( 12)-( 16).The FECVT model captures the cellularity present in the tissue microstructure, which in turn is used in the finite element analysis.In a conventional Voronoi model, the boundary cell areas are large because of semi-infinite edges.Hence, the areas of the virtual cells differ significantly with respect to the original images.The FECVT models of the selected collenchyma sections are simulated in ANSYS to determine the effective tissue stiffness.Figure 10 shows the nodal displacements of model tissue for 15% relative density with compressive strain along X and Y directions.The microstructure in X direction is less stiff than that in the Y direction.The microstructural anisotropy of the plant tissue originates from the cellular distribution.
We have compared the normalized stiffness obtained with the FEA of the FECVT model with the stiffness of a randomly generated Voronoi model, and the stiffness calculated from closed-form expressions available in literature [27].These formulas are obtained with an isotropic periodic cellular model having a hexagon as unit cell.On the other hand, the random Voronoi model is generated for a set of uniformly distributed points.Both the FECVT and randomly generated Voronoi models have the same number of cells.Figure 11 shows a regular hexagonal unit cell and the randomly generated Voronoi model.Both the randomly generated Voronoi and FECVT models are simulated for a set of relative densities between 5% and 30%.The normalized effective stiffness of both the models and the hexagonal unit cell are shown in Figure 12.The FEA of the random Voronoi model shows an average of 8% and 6% higher axial and shear modulus, respectively, compared with the closed-form solutions obtained for the unit cell shown in Figure 12.The results of the moduli are in agreement with those presented by Gibson et al. [26,27].For given relative density, the normalized effective elastic moduli of the different models in X and Y direction are shown in Figure 12.For relative density, 0.05 ≤ ρ * /ρ s ≤ 0.30, the FECVT model is 54% to 40% less stiff than the randomly generated Voronoi model and 31% to 35% less stiff than the periodic unit cell along the direction.Along the direction, the FECVT model is 19% to 15% less stiff than the random Voronoi model and 1% to 14% less stiff than the unit cell model for relative densities varying between 5% and 30%.To generate the Voronoi tessellation, we impose a uniform distribution of points.Hence, the randomly generated Voronoi model displays a uniform cell size, a factor that influences the stiffness properties, as shown in Figure 12.In the FECVT model, both shape and size of the cells vary significantly with respect to the random Voronoi model and the hexagonal cell model.Therefore, the shape and size of the cells affect the normalized stiffness, which varies with density.Similarly, the shear modulus of the FECVT model is lower than both the random Voronoi model and the unit cell model.The normalized effective stiffness of the FECVT model along the longitudinal direction is around 15% to 25% higher than that of the FECVT model along transverse direction for the range of density considered here.The variation reflects the stiffening effect of the cell shape, size, and the cellular distribution.It also depicts the anisotropic behaviour of the collenchyma tissue of the rhubarb petiole.
Since the FECVT model contains smaller cells than those of the actual tissue, the stiffness could be overestimated.In addition, the short cell walls modeled by the beam element may impose superfluous stiffness.Figure 13 depicts the effective Poisson's ratios for the FECVT, random Voronoi, and hexagonal unit cell models.The finite element analyses of both the FECVT and random Voronoi models exhibit marginal difference between the major and minor Poisson's ratio.Because of uniform cellular distribution, the Poisson's ratios of the random Voronoi model depend weakly on relative density.Hence, the microarchitecture of the tissue also affects the Poisson's ratio.
The effective stiffness of the cellular tissue can be obtained based on the computed wall stiffness in the previous hierarchy.For a relative density of 15%, the effective stiffness of the collenchyma tissue is 41.1 MPa and 11.02 MPa along X and Y direction, respectively.Instead of generic composition of dicotyledonous plant cell wall, the volume fractions and MFA of the model plant should be determined experimentally, and the computational stiffness would be more representative of and closer to the actual tissue stiffness.Moreover, the model is simplified, ignoring the presence of multilayers, if any.Despite the limitations, the approach is able to capture the stiffness at different hierarchical orders and the overall effective stiffness of the hierarchical cellular tissue.

Conclusions
Cellular solids are prevalent both in nature and in man-made engineering structures.Their microstructures govern the mechanical response at the macroscale.A micromechanical analysis that captures the realistic arrangement of the microstructure can thus provide insight into the overall macroscale behaviour of a cellular material.The goal of this work is to develop a multiscale hierarchical methodology to determine the stiffness of cellular plant tissue.In this paper, the stiffness properties at two orders of the tissue hierarchy have been determined by modeling structural details at their respective length scales.The upper and lower bounds of stiffness of the cell wall have been approximated via the classical micromechanical models of Voigt and Reuss.The microfibril angle and the stiffness bounds at the sub-cellular level have been used to compute the effective stiffness of the cell wall.We have shown that the MFA of the S 2 layer influences the overall stiffness of a multilayered secondary cell wall.The effect of MFA on the stiffness properties of single layered collenchyma cell wall of rhubarb petiole has also been captured at this order of hierarchy.The MFA governs the stiffness properties of both single and multilayered cell wall.At the cellular level, the tissue has been realistically modeled with the FECVT algorithm, which is capable of representing the actual cell distributions within the collenchyma tissue.The geometric model is representative of the cellular tissue and allows to numerically obtain the elastic properties of the cellular tissue.The heterogeneous cellular distribution obtained with the FECVT method can more realistically capture the compliant properties of plant tissue.

Figure 2 .
Figure 2. Hierarchical organization of Rhubarb petiole.The overall tissue stiffness depends on the stiffness properties of two hierarchical levels: n = 1, and 2. (©US Department of Energy Genome Programs [19].)

Figure 5 .
Figure 5. Simulated test for determining the effective stiffness properties.

Figure 6 .
Figure 6.Stiffness properties of different cell wall layers for varying MFA.

Figure 7 .
Figure 7. Poisson's ratio of different cell wall layers for varying MFA.

( a )Figure 8 .
Figure 8. Paraffin-embedded rhubarb petiole cross-section stained with TBO and imaged with light microscopy at 20× magnifications.Approximately 16 photos were stitched together to create this composite image.Ep-epithelium, Co-collenchyma cells, VB-vascular bundle, Pa-parenchyma cells are visible.Scale bar = 50 µm.(b) An enlarged view of a section of collenchyma tissue.

Figure 9 .
Figure 9. FECVT Model of the collenchyma tissue of rhubarb petiole with intermediate grey scale image.

Figure 10 .
Figure 10.Nodal displacement of FECVT model under uniaxial displacement boundary condition: (a) X component of nodal displacement; (b) Y component of nodal displacement.

( a )Figure 11 .
Figure 11.(a) Hexagonal unit cell; (b) Randomly generated Voronoi model for a set of 196 uniformly distributed points.

Figure 12 .
Figure 12.Normalized elastic modulus as a function of relative density.

Figure 13 .
Figure 13.Effective Poisson's ratios as a function of relative density.

Table 4 .
ANOVA for the FECVT models of collenchyma tissue.