Micromechanical model of biphasic biomaterials with internal adhesion: Application to nanocellulose hydrogel composites

Accepted Manuscript Micromechanical Model of Biphasic Biomaterials with Internal Adhesion: Application to Nanocellulose Hydrogel Composites Mauricio R. Bonilla, P. Lopez-Sanchez, M.J. Gidley, J.R. Stokes PII: DOI: Reference:

S1742-7061(15)30162-8 http://dx.doi.org/10.1016/j.actbio.2015.10.032 ACTBIO 3936

To appear in:

Acta Biomaterialia

Received Date: Revised Date: Accepted Date:

17 June 2015 12 October 2015 20 October 2015

Please cite this article as: Bonilla, M.R., Lopez-Sanchez, P., Gidley, M.J., Stokes, J.R., Micromechanical Model of Biphasic Biomaterials with Internal Adhesion: Application to Nanocellulose Hydrogel Composites, Acta Biomaterialia (2015), doi: http://dx.doi.org/10.1016/j.actbio.2015.10.032

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Micromechanical Model of Biphasic Biomaterials with Internal Adhesion: Application to Nanocellulose Hydrogel Composites Mauricio R. Bonilla1, P. Lopez-Sanchez2, M.J. Gidley2, J.R. Stokes1* 1

ARC Centre of Excellence in Plant Cell Walls, School of Chemical Engineering, The University of Queensland, Brisbane, 4072, Australia

2

ARC Centre of Excellence in Plant Cell Walls, Centre for Nutrition and Food Sciences,

Queensland Alliance for Agriculture and Food Innovation, The University of Queensland, Brisbane, 4072, Australia

*corresponding author: Associate Professor Jason Stokes School of Chemical Engineering The University of Queensland Brisbane, Australia Tel: +61 7 33654361 Email: [email protected]

1

Abstract The mechanical properties of hydrated biomaterials are non-recoverable upon unconfined compression if adhesion occurs between the structural components in the material upon fluid loss and apparent plastic behaviour. We explore these micromechanical phenomena by introducing an aggregation force and a critical yield pressure into the constitutive biphasic formulation for transversely isotropic tissues. The underlying hypothesis is that continual fluid pressure build-up during compression temporarily supresses aggregation. Once compression stops and the pressure falls below some critical value, internal aggregation occurs over a time scale comparable to the poroelastic time. We demonstrate this model by predicting the mechanical response of bacterial nanocellulose hydrogel composites, which are promising biomaterials and a structural mimetic for the plant cell wall. Cross-linking of cellulose by xyloglucan creates an extensional resistance and substantially increases the compressive modulus under large compression and densification. In comparison, incorporating non-crosslinking arabinoxylan into the hydrogel has little effect on its mechanics at the strain rates investigated. These results assist in elucidating the mechanical role of these polysaccharides in the complex plant cell wall structure. They also suggest xyloglucan is a suitable candidate to tailor the stiffness of nanocellulose hydrogels in biomaterial design, which includes modulating cell-adhesion in tissue engineering applications. The model and overall approach may be utilised to characterise and design a myriad of biomaterials and mammalian tissues, particularly those with a fibrillar structure.

2

1. Introduction The mechanical analysis of hydrogels through rheometric measurements offers a convenient avenue to characterise a myriad of biphasic materials of biological relevance such as plant cell walls, hydrated mammalian tissues, and biocompatible scaffolds. In particular, confined and unconfined stress/relaxation tests have become routine in the characterization of articular cartilage [1-3], cornea [4], white matter [5] and other soft tissues [6-8], as well as in the assessment of elastomeric hydrogels as potential replacements for damaged load bearing tissue [9, 10]. More recently, the introduction of a novel rheometric approach that combines compression-relaxation steps with small amplitude oscillatory shear (SAOS) provides a comprehensive picture of the material properties and permeability during densification of bacterial nanocellulose hydrogels [11]. Interpretation of stress/relaxation data requires a model capable of accounting for water redistribution within the sample, as this is expected to make a major contribution to the viscoelastic response. Most hydrated soft materials can be treated as a homogeneous mixture of a compressible solid matrix and a mobile interstitial fluid, and several theories for concomitant flowdeformation in porous media have been proposed to predict their viscoelastic behaviour. The biphasic model [1, 12] is one such theory and was originally proposed for the analysis of articular cartilage. It is straightforward to prove that the biphasic model corresponds to a simplification of Biot’s treatment for poroelastic materials [13] for the particular case of an incompressible fluid and an incompressible solid matrix. The theory has been extended to include sample anisotropy [14, 15], nonlinear elasticity [16, 17] and the intrinsic viscoelasticity of the solid phase, most notably by DiSilvestro et al. [2, 18]. While originally intended to mathematically describe articular cartilage, the fundamental assumptions of biphasic theory are fully applicable to hydrated soft tissues and hydrogels. Indeed, versions of the theory have been employed to describe the viscoelasticity of hydrogel-based contact 3

lenses [19], chondrocyte seeded agarose [20] and collagen [21] hydrogels for joint replacement, as well as PVA hydrogels [22, 23]. However, a simple approach for the analysis of anisotropic fibrillar materials such as nanocellulose hydrogels that are susceptible to fibre reorientation and aggregation during compression is not readily available. Biphasic fibrilreinforced models explicitly considering fibre reorientation and dispersion have been proposed for articular cartilage [3, 24], and the mechanics of fibre reinforced materials where fibres irreversibly reorient has been tackled in the context of arterial wall mechanics [25] and bacterial cellulose [26, 27]. However, current modelling approaches of reorientation and aggregation are too complex for routine analysis of experimental data due to the large number of mechanical parameters and the need for intricate global optimization schemes. In other cases, the effect of fluid pressure build-up and dissipation on the dynamics of fibre orientation is not included in the analysis. In this work, we propose a simple approach for the incorporation of reorientation/aggregation effects into the traditional biphasic theory, which is readily applicable to linear and nonlinear poroviscoelastic systems. A yield pressure is introduced in the stress-strain constitutive equations with a simple principle in mind: fibre-fibre adhesion forces are counteracted by continual pressure build up during compression, avoiding aggregation. However, once compression stops and fluid pressure falls below the yield pressure, aggregation occurs over a time scale that is expected to be similar to the poroelastic time. The proposed model is employed for the interpretation of unconfined compression of cellulose (C), cellulose – xyloglucan (CXG) and cellulose – arabinoxylan (CAX) hydrogels. We use cellulose hydrogels and its hemicellulosic composites as a simplified model system of the plant cell wall [28, 29]. Primary plant cell walls contain mainly cellulose, hemicelluloses and pectin. Xyloglucans are typically the main hemicellulose in primary cell walls of dicotyledonous plants, while heteroxylans (of which arabinoxylan is one example) 4

constitute the largest group of non-cellulosic polysaccharides in the primary cell walls of grasses, cereals and related commelinid species. These two hemicelluloses are known to interact in different ways with cellulose: xyloglucan will be embedded in cellulose fibres, on cellulose surface and crosslinking cellulose fibres, whereas arabinoxylan will be deposited on the surface on cellulose fibres [30, 31]. The cellulose/hemicellulose interactions are driven by hydrogen bonding and hydrophobic forces, with no evidence of significant contribution from electrostatic forces or the formation of covalent bonding [32, 33]. The characterization of these hydrogels using the current model provides insight on the specific contribution of each of these biopolymers to the hydrogel mechanics. Moreover, it shows that at different levels of hydration xyloglucan can significantly modify hydrogel stiffness with little variation in the permeability, which can be exploited in the development of cellulose – xyloglucan hydrogels for tissue engineering purposes [34].

2. Experimental procedures 2.1 Composite preparation Cellulosic composites were produced following the method described by Chanliaud et al. [35] and Mikkelsen et al. [36]. Briefly, the Gluconacetobacter xylinus frozen strain ATCC 53524 (Manassas, VA, USA) was cultivated in Hestrin and Schramm medium at pH 5. To produce cellulose/xyloglucan composites, a 1%(w/v) tamarind xyloglucan (Lot 100402, Megazyme International Ireland Ltd., County Wicklow, Ireland) solution was mixed with double concentrated Hestrin and Schramm medium (1:1) before inoculation, leading to a final xyloglucan concentration of 0.5% w/v. Similar methods and concentrations were used for the cellulose/arabinoxylan composites (medium viscosity wheat arabinoxylan, Lot 40302b, Megazyme International Ireland Ltd., County Wicklow, Ireland). Composites were cultivated statically at 30°C for 72 hours. After cultivation they were harvested and washed 6 times with 5

ice-cold water under agitation at 100 rpm, and stored in 0.02% NaN3 solution and kept at 4°C until further analysis. Further details can be found elsewhere [31]. The resulting samples were disks with average diameters of 40 mm and initial thicknesses: 1.3, 1.9 and 2.8 mm for CXG, CAX and C respectively. The initial polysaccharide concentration was ~1.4% for CXG and CAX, compared to ~0.7% for the cellulose-only hydrogel. The ratio of cellulose: xyloglucan was 1:0.4 in CXG, whereas the ratio of cellulose: arabinoxylan was 1:1.13 in CAX. SEM micrographs of the uncompressed hydrogels are reproduced in Figure 1 [31]. While XG cross-links cellulose fibrils, AX adsorbs onto the surface in discrete pockets. Such architectural arrangement of the hemicelluloses has a significant effect on the mechanical parameters of the different hydrogels, as will be shown below.

2.2 Compression-Relaxation test The mechanical testing methodology was described in a previous work for pure cellulose analysis [11], but is summarized here for completeness. Measurements were made on a rotational rheometer (HAAKE Mars III Rheometer, Thermo Fisher Scientific, Karlsruhe, Germany) at a constant temperature of 25°C controlled by a Peltier element. Parallel plates with a diameter of 60 mm were used. The upper and bottom plates were coated with fine Emery paper (P240/S85, 58 µm roughness). The samples were handled with tweezers and the composites were placed in the middle of the parallel plates with the help of a stencil. For each experiment the initial gap (distance between top and bottom plates) was adjusted according to the sample thickness. The behaviour under pressure is analysed here at two different compressive strain rates, 1 and 100 µm/s. Samples were compressed from their initial thickness in gap decrements of 100 6

µm, each followed by a relaxation period of 60 s. Compression continued to the narrowest possible gap before overshooting the normal forces allowed by the instrument. Between 6 and 9 replicates were measured. In this work, interpretation of the experimental data is performed after the samples are compressed to 1.6 mm, 1.1 mm and 1.3 mm for C, CAX and CXG, respectively, guaranteeing that all samples have an initial cellulose concentration of 0.0125 g/cm3. At this initial level of compression all samples relaxed completely and no residual normal pre-stress was registered. Moreover, no initial radial expansion was observed, suggesting absence of radial pre-stress. 3. Modelling of the micromechanical behaviour using biphasic theory According to Figure 1, the typical pore size of the hydrogels is somewhere between 0.1 and 1 µm, indicating that water holding is largely associated with capillary effects rather than physical or chemical adsorption. Hence, the equilibrium mechanical properties of the system (i.e., those measured when no fluid pressure gradient exists within the structure) are fully determined by the mechanical properties, concentration and architecture of the fibrillar network, with no interference from bulk properties of the fluid such as density or viscosity. However, if the time-scale of the experiment is comparable to the characteristic time of hydrodynamic dissipation, the time-dependent response of the material will be partially or completely dominated by fluid flow. In highly hydrated elastic materials, this characteristic time corresponds to the so-called poroelastic time tp [37, 38], which depends on both the mechanical properties and architecture of the fibre network and the fluid bulk properties. If t0 is the compression time and tr the measured relaxation time, then tp << t0 implies that hydrodynamic dissipation occurs very quickly and there will not be substantial viscoelastic effects associated to fluid flow. On the other hand, if tp >> tr, water transport is too slow and viscoelastic relaxation will not be significantly observed during the measurement period. For

7

a linear poroelastic material, tp ∼ µR2/kE¸ where µ is viscosity, k permeability, E the Young’s modulus and R the sample radius. A rough estimate of k is given by k ∼ 〈d2〉, where d is pore diameter. Our previous study determined that the compressive modulus for these cellulose hydrogels is between 10 4 and 105 Pa [11], which gives a tp ∼ 16 – 160s. The compression/relaxation experiments analysed here are all performed over time scales comparable to tp, which means that poroelasticity must be considered in addition the mechanics of the network.

We also highlight that poroelasticity is synonymous to

viscoelasticity arising from fluid flow through the porous structure. In sections 3.1 to 3.4, a new linear and nonlinear poroelastic model is introduced that incorporates fibre reorientation and aggregation. A more accurate formulation of the poroelastic time is introduced, but the magnitude of this characteristic time-scale remains the same as estimated above.

3.1 Linear biphasic model The composites are assumed to be transversely isotropic in the 1, 2-plane (Figure 2a). This is a reasonable consideration for the naocellulose hydrogels given that they are roughly produced in a layer-by-layer fashion, and have similarities to the configuration of the cellulosic network in plant cell walls [39, 40]. Transverse isotropy is also reasonable in the case of articular cartilage and other soft fibrilar tissues [4, 14, 41, 42]. In a previous work [11], an ad-hoc modification of the linear poroelastic theory developed by Cohen et al. [14] for the modelling of articular cartilage was successfully used to predict the response of pure nanocellulose hydrogels to unconfined compression/relaxation. According to this model, the normal stress σn(t) resulting from a ramp displacement in the z direction at constant strain rate ε
0 during t0 seconds, followed by a relaxation stage at constant strain is given by

8

σ n (t ) = E3ε
0t + E1

 1 ∞ exp(−α 2C kt / µ R 2 )  ∆ 3  − ∑ 2 2 2 n 11  C11k  8 i =1 α n α n ∆ 2 − ∆1 / (1 + ν 21 )  

µε
0 R 2

0 < t < t0

(1)

and

σ n (t ) = E3ε
0t0 + F *[1 − exp(−α12C11k (t − t0 ) / µ R 2 )]  ∞ exp(−α n2t / t g ) − exp[−α n2C11k (t − t0 ) / µ R 2 ]  + E1 ∆ 3 ∑  C11k α n2 α n2 ∆ 22 − ∆1 / (1 + ν 21 )   i=1 

µε
0 R 2

t > t0

(2)

where k is the permeability, µ the fluid viscosity, R the sample radius, E1 and ν12 the Young’s modulus and Poisson’s ratio in the plane of isotropy, respectively, and E3 and ν31 the out-of-plane Young’s modulus and Poisson’s ratio, respectively. αn corresponds to the roots of J1 ( x) − (1 −ν 312 E1 / E3 ) / (1 −ν 21 − 2ν 312 E1 / E3 ) xJ 0 ( x) = 0 , where J1 and J0 are Bessel functions of the first kind. Furthermore 2 E* ∆ 1 = 1 − ν 21 − ν 31

(3)

∆ 2 = (1 −ν 312 E * ) / (1 + ν 21 )

(4)

∆ 3 = (1 − 2ν 312 ) ∆ 2 / ∆ 1

(5)

E* = E1/E3. The aggregation modulus characterizing confined compression, C11, is given by:

C11 = E1 (1 −ν 312 E1 / E3 ) / [(1 + ν 21 )∆1 ]

(6)

The load F* < 0 (tensile) was superimposed to the normal stress to account for aggregation and reorientation of the fibres, at a rate that is assumed to be of the same order of magnitude 2 2 as that of water transport ∼ t p / α1 , where t p = µ R / C11k is the poroelastic time.

Bacterial cellulose fibres are pure 1-4 – β-glucans which do not have any surface charges. The fibres are not fundamentally adhesive in water, but can be forced into aggregation 9

through water depletion [31]. Neutron and X-ray scattering results [39] suggest that bacterial cellulose fibres in H2O display a core - shell structure. The core consists of cellulose crystallites surrounded by paracrystalline cellulose, while the shell is comprised by a significantly more porous layer of paracrystalline cellulose and tightly bounded water. Given that the shell thickness is about 4 nm and that its porosity may allow some degree of overlap between adjacent shells, their interactions are probably significant over a similar length-scale. Such interactions are likely the result of hydrogen bonding and dispersion forces.

Eq. (2) can be interpreted from the perspective of a pressure-dependent yield criterion: during compression fibre reorientation occurs, but fibre-fibre adhesion is counteracted by continual pressure build up, avoiding aggregation. However, once compression stops the fluid pressure dissipates at a rate ∼1/ tp, permitting the onset of aggregation over a time scale that is similar to tp. In the linear viscoelastic regime, this can be expressed through the constitutive formula

σ rr = C11ε rr + C12εθθ + C13ε zz − p σθθ = C12ε rr + C11εθθ + C13ε zz − p σ zz = C13ε rr + C13εθθ + C33 ( µ0 + µ1e

(7) −t /τ a

* zz

) ⊗ε − p

for p < py and ∂ p / ∂ t < 0 , and elsewhere by

σ rr = C11ε rr + C12εθθ + C13ε zz − p σθθ = C12ε rr + C11εθθ + C13ε zz − p σ zz = C13ε rr + C13εθθ + C33ε zz − p

(8)

Here, py is a yield pressure below which aggregation is significantly experienced during the relaxation stage. The ⊗ sign denotes the Stieltjes integral, t

g (t ) ⊗ f (t ) =

∫ g (t − τ ) ⊗ df (τ )

(9)

−∞

10

εii are the non-zero strain components, µ0 and µ1 are amplitude constants (µ0 + µ1 = 1), τa is an aggregation time constant and Cij’s are the stiffness matrix coefficients, related to the engineering material constants through C12 = E1 (ν 21 + ν 312 E * ) / [(1 + ν 21 ) ∆1 ] , C13 = E1ν 31 / ∆1 C33 = E3 [1 + 2ν 312 E * / ∆1 ] and eq. (6) for C11 [43].

εzz* = ε zz Θ( py − p) , where Θ(x) is the

Heaviside function taking the value of unity for x > 0 and zero elsewhere. Using the correspondence principle of viscoelasticity-elasticity, it can easily be proven [44, 45] that the solution to the unconfined compression problem for a ramp compressive strain is given by eqs. (1) and (2) when py is taken as equal or higher than the maximum pressure during the cycle, which for slow compression occurs at r = 0 and t = t0. The corresponding values of F* and τa are

F * = −C33µ1ε
0t0 ,

τ a = −µ R2 / α12C11k = t p / α12

(10)

*

Notice that the convolution in eq. (8) uses εzz (the instantaneous strain at p = py) as opposed to the strain history. It is reasonable to assume that fibre-fibre adhesion is a function of the average distance between fibres in contiguous layers along the axial direction, i.e., it is independent of the strain history prior to the time at which inter-fibre distance is small enough (below a critical distance dc) for adhesion to be significant. During compression, fibres in adjacent layers originally separated by an average distance above dc are drawn together (Figure 3a), but build-up of fluid pressure avoids adhesion even when the average distance between them falls below dc (Figure 3b). Pressure is dissipated in the relaxation stage but it is possible that initially the intrinsic adhesion force Fa is small compared to the hydraulic force Fh originated during fluid pressurization. Only when the hydraulic force is sufficiently low (below a certain yield pressure) does adhesion become observable on a macroscopic scale and an effective adhesion force F* manifests as a reduction in the axial 11

stiffness (Figure 3c). While simple, this approach provides a powerful tool for hydrogel characterization as will be shown below. Since the development of the biphasic model [1, 12], numerous modifications have been performed to introduce the effect of fibre reorientation during compression, particularly from the theory of fibre-reinforced composites [46]. In these works, either discrete spring elements representing the fibres are reoriented and deformed concomitantly with an isotropic elastic or viscoelastic matrix [47, 48] or both the matrix and springs are combined within a single continuous element [47, 49, 50]. While good agreement with experiment has been achieved for both large and small strains, these models are often difficult to implement and require complex optimization schemes to uniquely determine the model parameters, which can amount to ten or more. This has strongly limited their use for practical mechanical testing and development of hydrogels. The present approach allows straightforward implementation in standard programing software for linear systems or FEM solvers for nonlinear systems. The convenience of setting the effect of aggregation in the stiffness constant C33 is that, in the linear regime, it has no effect on the displacement field of the solid or the pressure distribution of the fluid for ideal unconfined compression. Consider the general case of

py < pmax , with pmax the maximum pressure during a cycle. The position r* at which p = py at time t during relaxation should be found prior to integration of the local normal stress over the sample surface. This implies inverting the pressure equation [1, 14] p=

 I 0 (r t p s ) − I 0 ( t p s ) C11 (C11 + C12 − 2C13 )    ε zz ( C11 − C12 )  C11 I 0 ( t p s ) − ( C11 − C12 ) I1 ( t p s ) / t p s 

(11)

from Laplace to time domain and solving for t at every r with p = py. I1 and I0 in eq. (11) are modified Bessel functions of the first kind and s is Laplace’s variable. A time domain solution of eq. (11) can easily be obtained using the residue theorem, 12

p (r , t ) =

(1 − 2υ31 ) E ∆ t ε
 ∞ J 0 (rα n / R ) / J 0 (α n ) − 1 exp(−α 2t / t )  n p 2 2 2 (1 + υ 21 ) 1 2 p 0 ∑ n =1 α n   ∆1 / (1 + υ 21 ) + ∆ 2α n    r  2   +0.64286 1 −      R   

(12)

t ≤ t0

(1 − 2υ 31 ) E ∆ t ε
× (1 + υ 21 ) 1 2 p 0  ∞ ( J 0 ( rα n / R ) / J 0 (α n ) − 1) ( exp( −α n2t / t p ) − exp( −α n2 (t − t 0 ) / t p ) )  ∑  α n2  ∆1 / (1 + υ 21 ) + ∆ 22α n2   n =1 

p(r , t ) =

(13) t > t0

Now, given that εrr =∂u / ∂r and εθθ = u / r one has, upon integration of σ zz in eqs. (7),

σ n (t ) =

2 R2

2µ C ε − 1 233 0 R

R

∫  − p + C (ε 13

rr

+ ε θθ ) + C 33ε zz  rdr

0 R

− ( t − t * ( r ) ) /τ a  ∫*  1 − e r (t )

(14)

  rdr 

where ε 0 = ε
0t0 . t*(r) corresponds to the time for which p = py at position r, while r*(t) corresponds to the position at which p = py at time t. Following the previous expression, transient radial gradients in material properties are expected to occur as portions of the sample near the edge experience a higher degree of aggregation during the relaxation stage. Static axial gradients in material properties are recognized to have an important impact on the mechanical response of soft materials [48, 51, 52] to external loads. Radial gradients and/or transient gradients in material properties have not been, to the best of our knowledge, investigated in detail. This will be part of future work. 3.2 Evolution of axial stiffness during compression/relaxation cycles

Figure 4 schematically depicts the evolution of the axial modulus along a series of compression/relaxation cycles as a function of (a) time and (b) compressive strain. Structural

13

changes in the microstructure of the composites produce an increase in the instantaneous axial modulus with compressive strain. On the other hand, when interpreting experimental data through the linear model above, one single “mean” axial modulus E3 is extracted (blue dots in Fig 5). However, during a single cycle the variation in axial modulus is presumably small if the change in compressive strain is small; hence taking a single value for E3 results in adequate fitting of the experimental data. During relaxation, reorientation and aggregation leads to a decrease in the instantaneous axial stiffness. Following eq. (2), the equilibrium f

modulus of the relaxed hydrogel E3 is given by

E3f = E3ε
0t0 +

F* = µ0 E3 ε
0t0

(15)

f

Since F* < 0, one has E3f < E3 . The equilibrium modulus E3 constitutes the starting stiffness for the next cycle, but the linear biphasic model provides a discrete increase in E3. If on the other hand the axial modulus is assumed to increase linearly with compressive strain the profile represented by the solid line and denoted as E3,nl in Figure 4 is obtained. In the next section, a piece-wise-linear and nonlinear scheme for the analysis of unconfined compression in hydrogels is introduced. 3.3 Piece-wise linear (PWL) poroviscoelasticity model

If nearly complete axial relaxation occurs up to relatively large nominal strains, it is reasonable to use the linear biphasic model in every compression/relaxation cycle, i.e., in a piece-wise linear fashion, to estimate the strain dependence of the system parameters. The normal stress is given by

σ n (t ) = ∆σ n (t ) + E3,f iε0,i

for ε 0,i−1 ≤ ε zz ≤ ε 0,i

(16)

14

f Where E3,i and ε0,i are equilibrium axial modulus and compressive strain at the end of cycle

i. ∆σn (t) is given by eqs. (1) and (2), with the axial strain (i.e. ε
0t in eq. (1) and ε
0t0 in eq. f (2)) replaced by ε zz − ε 0,i−1 . The equilibrium modulus E3, i is given by

E3,f i = ( µ0 E3 ) ε

(17) 0,i

where µ0 represents the fraction of the instantaneous axial modulus that does not vanish upon relaxation, following eq. (8). The piece-wise linear model will produce discrete average values of the fitting parameters, such as those represented by the symbols in Figure 4. 3.4 Non-linear biphasic model

In order to test the adequacy of the PWL approach, it is relevant to compare with a nonlinear model in the interpretation of experimental data. To introduce nonlinearity in the description of biphasic extracellular tissue, Wilson et al. [47] assumed that collagen fibrils in cartilage samples only resisted tension and were intrinsically viscoelastic, with hardening rates linearly related to the strain rate. An analogous approach has been used more recently in several works [3, 53] but requires testing at several strain rates for adequate determination of the viscoelastic parameters. Here we follow the approach by Li et al. [54], consisting of adding linear dependency of the moduli on the axial strain E1 = E10 + E1ε ε zz

E3 = E3,f i−1 + E3ε (ε zz − ε 0,i −1 )

(18)

for ε 0,i−1 ≤ ε zz ≤ ε 0,i

µ0 = µ00 + µ1ε εzz

(19)

(20)

15

The compressive strain is defined logarithmically as ε zz = ln ( h / h0 ) , with h and h0 the current and initial sample thicknesses, respectively. E1ε , E10 , µ00 and µ1ε are mechanical parameters to be determined. Notice that eq. (19) produces the profile in axial modulus denoted as E3,nl in Figure 4. The permeability is a function of the void ratio φ (pore volume over solid volume); the functional form below has proven effective for a number of hydrated tissues [47, 54, 55]

 φ − φ0  k = k 0 exp  M  1 + φ0  

(21)

where φ0 is the initial void ratio, taken as 0.98 as deduced from measurement of water content. We used the poroelastic module built into the COMSOL Multiphysics 4.4 to solve the equations of motion in the nonlinear case. An axisymmetric boundary is used at r = 0 (Figure 2). Contact with the platens is assumed frictionless with a roller boundary condition at z = 0. The discretization depicted in Figure 2b uses 420 quadratic elements, which were found to be adequate for the current problem. Solution with a finer mesh did not produce any significant improvement in the solution.

4. Results and Discussion 4.1 Cyclical slow compressions of nanocellulose hydrogels and composites

Figure 5a depicts the PWL and nonlinear fitting of compression/relaxation cycles of CAX and CXG at a compression rate of 1 µm/s, assuming py ≥ p0; similarly good fitting is obtained for C. The evolution of strain dependent parameters is depicted in Figure 6; in all scenarios,

ν12 > 0.9 was found to produce the most consistent fitting results. Each compression corresponds to a 100 µm decrease in gap (i.e., t0 = 100s). Direct measurement of the sample 16

diameter upon compression suggests that the equilibrium out-of-plane Poisson’s ratio is very small; consequently, υ31 was assumed to be zero for both the PWL and nonlinear models. For all hydrogels at strains below 30% it is noted that the high porosity of the samples produces a poroelastic time (which roughly corresponds to the relaxation time) that is small compared to t0, so that eq. (1) can be simplified as

σ n (t) = E3ε
0t +

ε
0 R2 8k

0 < t < t0

(22)

which means that E3 and k can be estimated from the slope m and intercept b of the linear section of the stress versus time curve just before t = t0. The difference between peak and relaxation stresses, ∆σ PR , yields

∆σ PR =

ε
0 R2 8k

− µ1E3ε
0t0

(23)

which allows direct extraction of µ1. If the relaxation time of the aggregation force τa is taken as tp/α12 [11], only two parameters are left to fit: E1 and υ21. At moderate and high strains one has tp ∼ t0 and only E3f can be directly determined by inspection of the residual stress. However, an informed initial guess of the remaining parameters can be made from the fitting performed at lower strains. The PWL model also produces a reasonable starting point for nonlinear fitting, reducing the risk of finding a physically irrelevant local minimum. The nonlinear prediction is slightly poorer than that from the PWL model, which is expected given that it contains a much lower number of fitting parameters to describe all cycles. Nevertheless, it performs remarkably well and further verifies the validity of the biphasic approach for the modelling of nanocellulose hydrogels and its composites. In all cases,

17

underprediction is observed at high strains, indicating that hardening occurs at a rate above that assumed by the linear hardening expressions in eqs. (18) - (19). Figures 6a and 6b depict the variation of the axial modulus with nominal compressive strain and cellulose density, respectively. The nonlinear model produces a profile in the axial modulus similar to that shown in Figure 4b; for simplicity, Figure 6a depicts only the maximum value at each ramp step. The fact that this maximum value increases nonlinearly with compressive strain is due to the use of finite deformations to define eqs. (19) and (20). Axial hardening is clearly significant, probably as a result of a larger number of fibres being bent, buckled or compressed to the point at which their contribution to axial stiffness shoots up nonlinearly. The agreement in the axial moduli extracted from the PWL and nonlinear models is remarkable, indicating that linear analysis of every cycle provides a good estimate of the evolution of the mechanical properties of this type of hydrogel. The good agreement is a consequence of the lack of residual radial stress upon relaxation, which means that only axial pre-stress (which can simply be added to the normal force) needs to be considered at the beginning of each cycle. When discrepancies due to cellulose concentration are ruled out (Figure 6b), it is clear that AX has essentially no effect on E3 at slow compression, while the presence of XG produces an increase of up to 30% in E3 with respect to the pure cellulosic hydrogel. The explanation for this is readily available from the micrographs in Figure 1: deposits of AX on the surface of cellulose fibrils do not provide any additional axial stiffness to the cellulose network, while the XG fibres crosslinking the cellulose fibrils are susceptible to compression, bending and buckling processes that boost axial stiffness. The XG fibres are however considerably weaker than the cellulose fibrils and their contribution manifests significantly only above 40% nominal strain. Interestingly, conventional theories of composite mechanics in soft tissues assume all compressive stiffness comes from a non-fibrilar isotropic matrix [46, 47, 18

54, 56]. Absence of such element in the cellulose hydrogels suggests that the contribution from fibre deformation is not necessarily negligible. There is a large difference in the magnitudes of the radial moduli (Figure 6c) for the different samples, with CXG displaying higher tensional stiffness than C and CAX. This apparently contradicts previous results from interpretation of uniaxial tension experiments that suggest xyloglucan significantly weakens tensional stiffness due to cellulose entanglement inhibition [57]. However, it must be considered that radial strains arising as a consequence of internal pressure gradients are infinitesimal, while uniaxial tension experiments involve network rearrangements on a much longer length-scale. Instead, the resulting radial stiffness in CXG is a combination of the strain applied on a collection of cellulose cross-links and that applied on XG fibres interconnecting cellulosic fibrils, as schematically shown in Figure 7. For all samples, the PWL and nonlinear models produce a radial modulus that increases with increasing compressive strain. The origin of in-plane hardening is likely to be due to the increasing volumetric density of cross-links with decreasing gap [11]. The same relationship between cross-links density and tensile hardening has been reported for polyethyleneglycol hydrogels [58]. Figures 6d depict the variation of permeability with nominal compressive strain. At slow compression rates there are no significant differences between the permeability of the pure cellulosic hydrogels and that of the composites. The structure of C is expectedly more open and allows easier water flow; however, its slightly higher permeability is in reality not enough to set its poroelastic time above that of CXG: on the contrary, tp for C increases monotonically from ∼50s to 200s, while tp for CXG increases monotonically from ∼50s to 150s. The reason for this is the higher radial modulus of CXG (tp ∼ 1/E1) compared to C, which allows for faster relaxation. In essence, these results suggest that XG not only 19

enhances compressive recovery, but also helps accelerate it. In contrast, AX does not affect the poroelastic time at slow compression rates because it does not affect significantly the architecture of the fibrillar network. The variation of the equilibrium axial modulus with nominal compressive strain following eq. (17) is depicted in Figure 6e. It is clear that the structure of all samples becomes progressively more elastic as the sample thickness decreases. SEM micrographs have previously suggested [31] that the structure of CXG and CAX hydrogels changes markedly when slowly compressed. In particular, fibril diameter and cross-linking density are increased, producing a structure that seems microscopically more robust. The equilibrium modulus of CXG exhibits a particularly rapid increase at high strains, suggesting that an additional mechanism is activated at small gaps. Under this level of confinement XG polymer chains might display intense repulsive interactions which manifest in hyperelastic hardening. For CXG hydrogels with lower water content (as e.g. in plant cell walls), XG is likely to improve elastic recovery upon compression. The increase in mechanical stiffness during densification is due to a combination of mechanical entanglement and adhesion. On the molecular scale, as there no charges on the cellulose fibres, both adhesion and friction between fibres is likely to be strongly determined by hydrogen bonding between the –OH groups on the surface of the fibres. The relative contribution of friction (tangential surface forces) and normal adhesion (perpendicular to the contact area) cannot be disentangled through the current model. However, ongoing research in our group is developing Atomic Force Microscopy based techniques that measure the adhesion and friction between fibres; these studies currently indicate that the energy required to separate a fibre from the network is similar to the adhesion energy arising from hydrogen bonding.

20

4.2 Fast compression

Unconfined compression/relaxation tests at a rate of 100µm/s were performed on the different composites; the results obtained from interpreting the experimental data with the PWL model for CAX are depicted in Figure 5b (left). Similar degree of agreement between experiment and model was obtained for C and CXG. In all cases the yield pressure py was below the maximum pressure: py = 185 Pa for C, py = 202 Pa for CXG and py = 235 Pa for CAX. Setting a finite py below the maximum pressure produces better curve fitting, implying that aggregation occurs only after part of the hydraulic pressure has been dissipated. While py is potentially a function of strain, for simplicity we have used a single, strain-independent py for our model fitting. In all cases py is above the maximum pressure achieved during slow compression for the hydrogels, explaining why the use of a critical pressure is not required in the slow compression scenario. Figure 5b (right) depicts the effect of making py larger than the maximum pressure in one of the compression/relaxation cycles during of CAX. Clearly, assuming aggregation to occur from the onset of relaxation does not provide adequate interpretation of the experimental results. The values of py between the different composites do not differ significantly, suggesting that aggregation is mostly driven by direct cellulosecellulose interactions. Most models of cellulose hydrogel composite systems like plant cell walls disregard the mechanical effect of direct cellulose interactions, assuming cross-linking polymers such as XG are entirely responsible for stress transfer between cellulose fibres [5962]. These models have been challenged [27, 57], with alternative models proposing that XG and AX act mostly as modulators of the cellulose network architecture, while direct entanglement and adhesion between cellulose fibres constitute the main stress transfer mechanism. While the current results do not directly support the latter model, it shows that at the investigated compositions the hydraulic pressure necessary to inhibit the onset of adhesion is almost independent XG and AX concentrations. 21

Fitting of slow and fast compression data produces different values of the mechanical parameters, though within the same order of magnitude and with the similar variation with compressive strain. In particular, higher values of the instantaneous axial modulus E3 are attained during fast compression (not shown) compared to slow compression, which means that the systems display intrinsic viscoelasticity in the axial direction. However, we have refrained here from proposing any particular form for the intrinsic viscoelasticity given that proper validation would require experimental data at a comprehensive number of compression speeds. This is straightforward to do by, for example, replacing the coefficients in eqs. (7) and (8) by appropriate time decaying functions and the coefficient-strain multiplications by convolutions, in the spirit of Hoang et al. [44, 45] Moreover, previous work on cellulose composites has shown that significant structural differences exist between fast and slowly compressed samples [31] and consequently, the material properties are expected to differ. If py is not introduced an artificially high value for E3 is obtained, because it is not possible to reconcile the fast initial increase in normal force during compression and its slow decay at the end of the relaxation stage with a unique value of permeability. In a way, py allows the system to “regulate” where the slow decay begins, while the decay rate is controlled through the aggregation time constant τa. There is general consistency in the relative order of E1 and E3 for the different hydrogels: CXG is radially and axially stiffer than CAX and C for both slow and fast compression. However, such consistency is absent regarding the permeability (Figures 8). While the permeability of the three samples extracted from slow compression experiments is essentially the same, fitting of fast compressions shows that CAX has lower permeability than the other samples. The discrepancy suggests that the large pressure gradients generated during fast

22

compression causes the AX to effectively increase the fluid viscosity. This is consistent with the finding that the AX associated with cellulose fibrils retains segmental flexibility as shown by NMR spectroscopy [32]; under confined conditions this is expected to generate high local viscosity. As a matter of fact, CAX displays considerably higher peak forces than C when compressed at high strain rate, despite having at the starting point the same cellulose concentration.

23

5. Concluding Remarks

In this work, we propose a modified biphasic theory to interpret the micromechanical behaviour of nanocellulose hydrogel composites under unconfined compression. The main hypothesis is that fluid pressure dissipation controls the rate of aggregation, which is responsible for producing equilibrium axial moduli below that expected from the dynamics of the compression stage. Once compression stops and the fluid pressure begins to dissipate, aggregation is assumed to occur only after the hydraulic force counteracting the adhesive interaction between fibres falls below a certain pressure threshold. Such an effect is more obvious during fast compressions, where fluid pressure increases more significantly. In the absence of residual radial stress, it is shown that a piece-wise linear approach suffices to predict multiple unconfined compression/relaxation cycles, avoiding the need for complex nonlinear models and optimization schemes for which a unique set of mechanical parameters is generally not guaranteed. Indeed, direct observation of the compression/relaxation curves for compression times larger than the characteristic relaxation time provides an easy route for estimating the material parameters. The model introduced here is relevant for the analysis of hydrated mammalian tissues such as articular cartilage or nucleus pulposus, given their structural similarity with the hydrogels in this study. Moreover, the high crystallinity, water holding capacity and high tensile strength of nanocellulose hydrogels and associated composites make them viable materials for bone and cartilage scaffolding as well as muscle and blood vessel reconstruction. Availability of suitable models for mechanical analysis of hydrogel-like tissues and composites is consequently of paramount importance for rational assessment of new biomaterials for tissue engineering applications; we have contributed to this endeavour by providing a model of relatively straightforward implementation.

24

While cellulose is known to provide extensional stiffness, it does not provide significant compressive stress recovery. AX does not have any direct mechanical effect when present in the hydrogel, consistent with physical adsorption onto the surface of the cellulose fibrils. However, the residual segmental flexibility of AX means that it is possible that at much higher concentration it produces an interconnected isotropic matrix that can receive and transmit load to the fibrilar structure, producing a fibril reinforced composite in the traditional sense. Since XG chains bridge cellulose fibrils, this provides extensional resistance to the tensional strains generated by internal pressure gradients, which are generally small compared to those experienced during uniaxial extensional testing. Since there is no amorphous matrix sustaining the load upon repeated compressions, the compression/relaxation test does not produce radial stresses that are strong enough to probe all cellulosic crosslinks. Therefore, the effect of XG crosslinks has a relevant role in increasing the radial modulus compared to the purely cellulosic hydrogel in compression/relaxation experiments. XG also produces an increase of the compressive modulus at large strains, which translates into an overall increase in the resistance to compression of the hydrogels. XG has been previously used to enhance cell adhesion onto cellulose – based scaffolds and its potential for mechanical modification of nanocellulose hydrogels for biomimetic purposes has been recognized. However, the effect XG incorporation has on the response to compressional forces had not been addressed before on a quantitative fashion. While this work does not attempt to be comprehensive on hydrogel tailoring through XG addition, it shows that its presence can significantly improve the compressive stiffness of cellulose scaffolds, particularly at large strains with almost no variation in hydrogel permeability compared to purely nanocellulose hydrogels. This means XG could lead to structures that can better withstand very large strains while avoiding diffusional slow-down of nutrient 25

exchange. Future work should focus on achieving prescribed mechanical properties through xyloglucan incorporation.

ACKNOWLEDGMENTS The authors would like to thank Dr. Deirdre Mikkelsen for

helpful advice during the production of bacterial cellulose composites used for mechanical testing. The authors thank Grace Dolan for helpful discussions and her input concerning AFM measurements on adhesion and friction between fibres. The research has been funded by the Australian Research Council Centre of Excellence on Plant Cell Walls.

26

Figure Captions Figure. 1 Scanning electron micrographs of a) cellulose (C), b) cellulose/xyloglucan (CXG)

and c) cellulose/arabinoxylan (CAX) hydrogels. Xyloglucan crosslinks are present in CXG composite. Arabinoxylan forms discrete aggregates onto the cellulose fibres. From LopezSanchez et al. [31]. Figure 2. (a) Geometry of the hydrogel sample, where 1-2 is the plane of transverse isotropy.

Axisymmety allows reducing the problem to two dimensions. (b) Discretization of an axisymmetric slice of the hydrogel disk using quadratic elements in Comsol multiphysics 4.4. Figure 3. Fibre-fibre adhesion force Fa is a function of the average distance between fibres d

and is only significant if d is below a certain distance threshold dc. During compression, fibres in contiguous layers are drawn together, but build-up of fluid pressure avoids adhesion (a). During relaxation, pressure is dissipated and the average inter-fibre distance decreases below dc; however, it is possible that the intrinsic adhesion force Fa is small compared to the hydraulic force Fh, which is itself a function of pressure. In this case, the macroscopic effect of the adhesion force F* ∝ Fa − Fh is negligible (b). Only for Fh(p) ≤ Fh(py) does F* becomes significant enough to be macroscopically registered (c). Figure 4. Evolution of the axial modulus along a series of compression/relaxation cycles as a

function of (a) time and (b) compressive strain. Microstructural changes in the hydrogels leads to increasing instantaneous axial modulus with increasing compressive strain. Interpreting experimental data through the linear poroelastic model produces an “average” axial modulus E3 (blue dots). During relaxation, reorientation and aggregation leads to a decrease in axial stiffness. The equilibrium modulus of the relaxed hydrogel E 3f is given by

27

eq. (15). E 3f constitutes the starting stiffness for the next cycle, but the linear model provides a discrete increase in E3. If the axial modulus is assumed to increase linearly with compressive strain (E3,nl) the profile represented by the solid line is obtained. Figure 5. (a) PWL and nonlinear fitting of compression/relaxation cycles of CAX (left) and

CXG (right) at a compression speed of -1µm/s, assuming py > p(t0, 0). While both models show good agreement with experiments, some underestimation is observed at high strains in the normal stress predicted by the nonlinear model, as a consequence of the assumed shape of the nonlinearities in eqs. (18) -(20). Similar level of agreement between model and experiment was obtained for C at this compression speed. (b) PWL fitting of compression/relaxation cycles of CAX (left) at a compression speed of -100 µm/s, assuming constant py throughout the experiment. The single-cycle zoom on the right depicts the effect of making py larger than the maximum pressure in one of the cycles in CAX. Clearly, assuming aggregation to occur from the onset of relaxation does not adequately fit the experimental results. Figure 6. Mechanical parameters obtained from interpretation of slow unconfined

compression experiments through the PWL and nonlinear models. Figure 7. CXG displays higher tensional stiffness than C and CAX. This is due to cross-

linking of cellulose fibrils by XG fibres: internal pressure gradients tend to separate cellulose fibrils, which induces tensional strain on the XG fibres. This produces additional radial stiffness that ultimately translates into a higher radial modulus for CXG. Figure 8. Permeability values obtained from interpretation of fast unconfined compression

experiments through the PWL model.

28

References

[1] Armstrong CG, Lai WM, Mow VC. An Analysis of the Unconfined Compression of Articular Cartilage. Journal of Biomechanical Engineering 1984;106:165-73. [2] DiSilvestro MR, Zhu Q, Suh J-KF. Biphasic Poroviscoelastic Simulation of the Unconfined Compression of Articular Cartilage: II—Effect of Variable Strain Rates. Journal of Biomechanical Engineering 2000;123:198-200. [3] Seifzadeh A, Wang J, Oguamanam DCD, Papini M. A Nonlinear Biphasic FiberReinforced Porohyperviscoelastic Model of Articular Cartilage Incorporating Fiber Reorientation and Dispersion. Journal of Biomechanical Engineering 2011;133:081004-. [4] Hatami-Marbini H, Etebu E. Rate dependent biomechanical properties of corneal stroma in unconfined compression. Biorheology 2013;50:133-47. [5] Cheng S, Bilston LE. Unconfined compression of white matter. Journal of Biomechanics 2007;40:117-24. [6] Wang P, Yang L, Hsieh A. Nucleus Pulposus Cell Response to Confined and Unconfined Compression Implicates Mechanoregulation by Fluid Shear Stress. Ann Biomed Eng 2011;39:1101-11. [7] Cloyd JM, Malhotra NR, Weng L, Chen W, Mauck RL, Elliott DM. Material properties in unconfined compression of human nucleus pulposus, injectable hyaluronic acid-based hydrogels and tissue engineering scaffolds. European Spine Journal 2007;16:1892-8. [8] Miller K. Method of testing very soft biological tissues in compression. Journal of Biomechanics 2005;38:153-8. [9] Liu K, Ovaert TC. Poro-viscoelastic constitutive modeling of unconfined creep of hydrogels using finite element analysis with integrated optimization method. Journal of the Mechanical Behavior of Biomedical Materials 2011;4:440-50. [10] Lee S-Y, Pereira BP, Yusof N, Selvaratnam L, Yu Z, Abbas AA, et al. Unconfined compression properties of a porous poly(vinyl alcohol)–chitosan-based hydrogel after hydration. Acta Biomaterialia 2009;5:1919-25. [11] Lopez-Sanchez P, Rincon M, Wang D, Brulhart S, Stokes JR, Gidley MJ. Micromechanics and Poroelasticity of Hydrated Cellulose Networks. Biomacromolecules 2014;15:2274-84. [12] Mow VC, Kuei SC, Lai WM, Armstrong CG. Biphasic Creep and Stress Relaxation of Articular Cartilage in Compression: Theory and Experiments. Journal of Biomechanical Engineering 1980;102:73-84. [13] Biot MA. General theory of three dimensional consolidation. Journal of Applied Physics 1941;12:155–64. [14] Cohen B, Lai WM, Mow VC. A Transversely Isotropic Biphasic Model for Unconfined Compression of Growth Plate and Chondroepiphysis. Journal of Biomechanical Engineering 1998;120:491-6. [15] Soltz MA, Ateshian GA. A Conewise Linear Elasticity Mixture Model for the Analysis of Tension-Compression Nonlinearity in Articular Cartilage. Journal of Biomechanical Engineering 2000;122:576-86. [16] Garcia JJ, Altiero NJ, Haut RC. Estimation of in Situ Elastic Properties of Biphasic Cartilage Based on a Transversely Isotropic Hypo-Elastic Model. Journal of Biomechanical Engineering 1999;122:1-8.

29

[17] Fortin M, Soulhat J, Hunziker EB, Buschmann MD, Shirazi-Adl A. Unconfined Compression of Articular Cartilage: Nonlinear Behavior and Comparison With a FibrilReinforced Biphasic Model. Journal of Biomechanical Engineering 1999;122:189-95. [18] DiSilvestro MR, Wong M, Jurvelin JS, Suh J-KF, Zhu Q. Biphasic Poroviscoelastic Simulation of the Unconfined Compression of Articular Cartilage: I—Simultaneous Prediction of Reaction Force and Lateral Displacement. Journal of Biomechanical Engineering 2000;123:191-7. [19] Chen X, Dunn AC, Sawyer WG, Sarntinoranont M. A Biphasic Model for MicroIndentation of a Hydrogel-Based Contact Lens. Journal of Biomechanical Engineering 2006;129:156-63 [20] Kim E, Guilak F, Haider MA. The dynamic mechanical environment of the chondrocyte: A biphasic finite element model of cell-matrix interactions under cyclic compressive loading. Journal of Biomechanical Engineering 2008;130:061009-. [21] Galie P, Spilker R, Stegemann J. A Linear, Biphasic Model Incorporating a Brinkman Term to Describe the Mechanics of Cell-Seeded Collagen Hydrogels. Ann Biomed Eng 2011;39:2767-79. [22] Goldsmith AAJ, Clift SE. Investigation Into the Biphasic Properties of a Hydrogel for Use in a Cushion Form Replacement Joint. Journal of Biomechanical Engineering 1998;120:362-9. [23] Blum MM, Ovaert TC. Experimental and numerical tribological studies of a boundary lubricant functionalized poro-viscoelastic PVA hydrogel in normal contact and sliding. Journal of the Mechanical Behavior of Biomedical Materials 2012;14:248-58. [24] Quinn TM, Morel V. Microstructural Modeling of Collagen Network Mechanics and Interactions with the Proteoglycan Gel in Articular Cartilage. Biomech Model Mechanobiol 2007;6:73-82. [25] Holzapfel G, Gasser T, Ogden R. A New Constitutive Framework for Arterial Wall Mechanics and a Comparative Study of Material Models. Journal of Elasticity 2000;61:1-48. [26] Gindl W, Keckes J. Tensile properties of cellulose acetate butyrate composites reinforced with bacterial cellulose. Composites Science and Technology 2004;64:2407-13. [27] Astley OM, Chanliaud E, Donald AM, Gidley MJ. Tensile deformation of bacterial cellulose composites. International Journal of Biological Macromolecules 2003;32:28-35. [28] Chanliaud E, De Silva J, Strongitharm B, Jeronimidis G, Gidley MJ. Mechanical effects of plant cell wall enzymes on cellulose/xyloglucan composites. Plant J 2004;38:27-37. [29] McKenna BA, Mikkelsen D, Wehr JB, Gidley MJ, Menzies NW. Mechanical and structural properties of native and alkali-treated bacterial cellulose produced by Gluconacetobacter xylinus strain ATCC 53524. Cellulose 2009;16:1047-55. [30] Pauly M, Albersheim P, Darvill A, York WS. Molecular domains of the cellulose/xyloglucan network in the cell walls of higher plants. The Plant Journal 1999;20:629-39. [31] Lopez-Sanchez P, Cersosimo J, Wang D, Flanagan B, Stokes JR, Gidley MJ. Poroelastic Mechanical Effects of Hemicelluloses on Cellulosic Hydrogels under Compression. PLoS ONE 2015;10:e0122132. [32] Mikkelsen D, Flanagan BM, Wilson SM, Bacic A, Gidley MJ. Interactions of Arabinoxylan and (1,3)(1,4)-β-Glucan with Cellulose Networks. Biomacromolecules 2015;16:1232-9. [33] Lopez M, Bizot H, Chambat G, Marais M-F, Zykwinska A, Ralet M-C, et al. Enthalpic Studies of Xyloglucan−Cellulose Interactions. Biomacromolecules 2010;11:1417-28. [34] Bodin A, Ahrenstedt L, Fink H, Brumer H, Risberg B, Gatenholm P. Modification of Nanocellulose with a Xyloglucan–RGD Conjugate Enhances Adhesion and Proliferation of

30

Endothelial Cells: Implications for Tissue Engineering. Biomacromolecules 2007;8:3697704. [35] Chanliaud E, Gidley MJ. In vitro synthesis and properties of pectin/Acetobacter xylinus cellulose composites. Plant J 1999;20:25-35. [36] Mikkelsen D, Gidley M. Formation of Cellulose-Based Composites with Hemicelluloses and Pectins Using Gluconacetobacter Fermentation. In: Popper ZA, editor. The Plant Cell Wall: Humana Press; 2011. p. 197-208. [37] Moeendarbary E, Valon L, Fritzsche M, Harris AR, Moulding DA, Thrasher AJ, et al. The cytoplasm of living cells behaves as a poroelastic material. Nat Mater 2013;12:253-61. [38] Skotheim JM, Mahadevan L. Physical Limits and Design Principles for Plant and Fungal Movements. Science 2005;308:1308-10. [39] Cave ID. The anisotropic elasticity of the plant cell wall. Wood Science and Technology 1968;2:268-78. [40] Baskin TI. Anisotropic Expansion of the Plant Cell Wall. Annual Review of Cell and Developmental Biology 2005;21:203-22. [41] Yin L, Elliott DM. A biphasic and transversely isotropic mechanical model for tendon:: application to mouse tail fascicles in uniaxial tension. Journal of Biomechanics 2004;37:90716. [42] Federico S, Grillo A, La Rosa G, Giaquinta G, Herzog W. A transversely isotropic, transversely homogeneous microstructural-statistical model of articular cartilage. Journal of Biomechanics 2005;38:2008-18. [43] Jones RM. Mechanics of Composite Materials. Washington, DC.: Scripta; 1975. [44] Hoang SK, Abousleiman YN. Poroviscoelasticity of transversely isotropic cylinders under laboratory loading conditions. Mechanics Research Communications 2010;37:298-306. [45] Hoang SK, Abousleiman YN. Poroviscoelastic Two-Dimensional Anisotropic Solution with Application to Articular Cartilage Testing. Journal of Engineering Mechanics 2009;135:367-74. [46] Spencer AJM. Continuum Theory of the Mechanics of Fibre-Reinforced Composites. Wien: Springer,; 1984. [47] Wilson W, van Donkelaar CC, van Rietbergen B, Huiskes R. A fibril-reinforced poroviscoelastic swelling model for articular cartilage. Journal of Biomechanics 2005;38:1195-204. [48] Li LP, Cheung JTM, Herzog W. Three-dimensional fibril-reinforced finite element model of articular cartilage. Med Biol Eng Comput 2009;47:607-15. [49] Julkunen P, Kiviranta P, Wilson W, Jurvelin JS, Korhonen RK. Characterization of articular cartilage by combining microscopic analysis with a fibril-reinforced finite-element model. Journal of Biomechanics 2007;40:1862-70. [50] Gasser TC, Ogden RW, Holzapfel GA. Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. Journal of The Royal Society Interface 2006;3:15-35. [51] Armstrong CG, Mow VC. Variations in the intrinsic mechanical properties of human articular cartilage with age, degeneration, and water content1982. [52] Calve S, Simon H-G. Biochemical and mechanical environment cooperatively regulate skeletal muscle regeneration. The FASEB Journal 2012;26:2538-45. [53] García JJ, Cortés DH. A nonlinear biphasic viscohyperelastic model for articular cartilage. Journal of Biomechanics 2006;39:2991-8. [54] Li LP, Soulhat J, Buschmann MD, Shirazi-Adl A. Nonlinear analysis of cartilage in unconfined ramp compression using a fibril reinforced poroelastic model. Clinical Biomechanics 1999;14:673-82. [55] Argoubi M, ShiraziAdl A. Poroelastic creep response analysis of a lumbar motion segment in compression. Journal of Biomechanics 1996;29:1331-9. 31

[56] Soulhat J, Buschmann MD, Shirazi-Adl A. A Fibril-Network-Reinforced Biphasic Model of Cartilage in Unconfined Compression. Journal of Biomechanical Engineering 1999;121:340-7. [57] Chanliaud E, Burrows KM, Jeronimidis G, Gidley MJ. Mechanical properties of primary plant cell wall analogues. Planta 2002;215:989-96. [58] Nguyen QT, Hwang Y, Chen AC, Varghese S, Sah RL. Cartilage-like mechanical properties of poly (ethylene glycol)-diacrylate hydrogels. Biomaterials 2012;33:6682-90. [59] Kha H, Tuble SC, Kalyanasundaram S, Williamson RE. WallGen, Software to Construct Layered Cellulose-Hemicellulose Networks and Predict Their Small Deformation Mechanics. Plant Physiology 2010;152:774-86. [60] Yi H, Puri VM. Contributions of the mechanical properties of major structural polysaccharides to the stiffness of a cell wall network model. American Journal of Botany 2014;101:244-54. [61] Yi H, Puri VM. Architecture-Based Multiscale Computational Modeling of Plant Cell Wall Mechanics to Examine the Hydrogen-Bonding Hypothesis of the Cell Wall Network Structure Model. Plant Physiology 2012;160:1281-92. [62] Zykwinska AW, Ralet M-CJ, Garnier CD, Thibault J-FJ. Evidence for In Vitro Binding of Pectin Side Chains to Cellulose. Plant Physiology 2005;139:397-407.

32

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Graphical abstract

33