Equivalent mechanical properties of biological membranes from lattice homogenization

J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S

4 (2011) 1833–1845

Available online at www.sciencedirect.com

journal homepage: www.elsevier.com/locate/jmbbm

Research paper

Equivalent mechanical properties of biological membranes from lattice homogenization M. Assidi 1 , F. Dos Reis, J.-F. Ganghoffer ∗ Laboratoire d’Énergétique et de Mécanique Théorique et Appliquée, École Nationale Supérieure d’Électricité et de Mécanique, UMR 7563, ENSEM-INPL, 2, Avenue de la Forêt de Haye, BP 160, 54054 Vandoeuvre CEDEX, France

A R T I C L E

I N F O

A B S T R A C T

Article history:

The goal of this manuscript is to set up a novel methodology for the calculation of the

Received 22 March 2011

effective mechanical properties of biological membranes viewed as repetitive networks of

Received in revised form

elastic filaments, based on the discrete asymptotic homogenization method. We will show

20 May 2011

that for some lattice configurations, flexional effects due to internal structure mechanisms

Accepted 28 May 2011

at the unit cell scale lead to additional flexional effects at the continuum scale, accounted

Published online 12 June 2011

for by an internal length associated to a micropolar behavior. Thereby, a systematic methodology is established, allowing the prediction of the overall mechanical properties

Keywords:

of biological membranes for a given network topology, as closed form expressions of the

Biological membranes

geometrical and mechanical micro-parameters. The peptidoglycan and the erythrocyte

Asymptotic homogenization

have been analyzed using this methodology, and their effective moduli are calculated and

Equivalent properties

recorded versus the geometrical and mechanical lattice parameters. A classification of

Nonlinear response

lattices with respect to the choice of the equivalent continuum model is proposed: The

Micropolar continuum

Cauchy continuum and a micropolar continuum are adopted as two possible effective medium, for a given beam model. The relative ratio of the characteristic length of the micropolar continuum to the unit cell size determines the relevant choice of the equivalent medium. In most cases, the Cauchy continuum is sufficient to model membranes in most of their configurations. The peptidoglycan network may exhibit a re-entrant hexagonal lattice, for which micropolar effects become important. This is attested by the characteristic length becoming larger than the beam length for such configurations. The homogenized moduli give accurate results for both membranes, as revealed by comparison with experimental measurements or simulation results from the literature at the network scale. A first insight into the nonlinear mechanical behavior of the hexagonal and triangular networks is lastly investigated using a perturbative method. c 2011 Elsevier Ltd. All rights reserved. ⃝

1.

Introduction

The membrane of biological cells is made of the assembly of filaments which are linked together as part of a network

or are associated with the cell membrane to build a twodimensional thin sheet. Two-dimensional biological networks may be wrapped around a cell as its wall or attached to its plasma or nuclear membrane. Structural elements

∗ Corresponding address: ENSEM-INPL, 2 Avenue de la Forêt de Haye, 54500 Vandoeuvre-Les-Nancy, France. Tel.: +33 3 83 59 57 24; fax: 33 3 83 59 55 51. E-mail addresses: [email protected] (M. Assidi), [email protected] (F. Dos Reis), [email protected] (J.-F. Ganghoffer). 1 Tel.: +33 621345616. c 2011 Elsevier Ltd. All rights reserved. 1751-6161/$ - see front matter ⃝ doi:10.1016/j.jmbbm.2011.05.040

1834

J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S

of biological cells are soft and responsible for the large deformability and easy motion of the cell, contrary to most of the engineered man-made thin structural materials. The mechanics of biological membranes is clearly related to the network architecture and to the elasticity of the building entities of the network, which are known to obey entropic elasticity. The main originality advocated in this work (compared to the literature works referenced before) is the explicit derivation of the full stiffness matrix of biological membranes, viewed as planar repetitive networks of connected threads, as closed form expressions of the geometrical and mechanical micro-parameters of the underlying network. Moreover, we shall justify the use of continuum models with an enriched kinematics such as micropolar continua according to length scale considerations. As a third novel aspect, the construction of the effective mechanical response in the large strain regime using the same homogenization technique will be done. The development of predictive nanomechanical models aiming at understanding the impact of the network architecture and mechanical properties on the continuum scale is important, as the experimental determination of the mechanical properties of biological membranes is delicate, and the membrane anisotropy and large deformability have to be accounted for, Boey et al. (1998), Discher et al. (1998). As stated in Lim et al. (2006), mechanical models for cells are derived using either the micro/nanostructural approach or the continuum one. Although providing less insight into detailed molecular mechanical events and biochemical couplings, the continuum approach is easier and more straightforward to use in computing the mechanical properties of the cell and its response under biomechanical loading. Moreover, the established continuum mechanical model can provide details on the distribution of stress and strains induced in the cell and can be integrated in finite element simulations at the scale of the whole cell. However, the identification of the continuum behavior of a membrane is challenging, as it may be highly anisotropic due to unequal chain length and properties of the threads within the molecular network may vary; furthermore, biological membranes are prone to large distensions and one should ideally consider nonlinear effects. Hence, micromechanical approaches are needed in order to bridge the scales and to provide a constitutive law at a continuum scale, whereby the equivalent continuum properties are related to both the geometrical and mechanical nanostructural parameters of the network. The derivation of the equivalent mechanical properties of cellular biological structures is also interesting in order to understand the somewhat peculiar observed behavior (anisotropy, negative Poisson’s ratio, (Boal et al., 1993)) and to possibly evaluate the load bearing capacity of the membrane architecture. Especially, closed form expressions of those effective properties would allow relating the mesoscopic to the maroscopic level, to understand the nanoscale origin of the mechanical behavior of the membrane wall, and to assess the effect of the membrane topology (comparison of different membrane architectures will be possible). We shall employ the so-called discrete asymptotic homogenization technique as in Caillerie et al. (2006), which is perfectly suited to the discrete architecture of the membrane at

4 (2011) 1833–1845

the nano-level. Two types of equivalent continuum shall be considered, a classical Cauchy continuum and a micropolar medium, according to the value of a characteristic micropolar length. This last aspect constitutes the main and novel thrust of this contribution, especially when considering biological membranes.

2.

Impact of microstructural irregularity

Homogenization techniques for discrete media have been extensively used in the last decade, but they have a significant limitation in that they do not account for natural variations in the lattice topology, which are observed for most biological materials. Most models of 2D cellular structures are based on idealized unit cells intended to describe the micro-structural features of an average cell supposed to be representative of the real underlying structure. Those approaches do however not account for the complex and rather diverse mechanisms leading to membrane rearrangements usually referred to as remodeling; those mechanisms involving a complex machinery of proteins can be broadly classified as fusion or fission, including exocytosis and endocytosis, budding and fusion of transport carriers, relaxation of the elastic energy, as listed in the recent review paper (Kozlov et al., 2010). The network topology may also vary as abnormal RBC skeletons have been reported, Hansen et al. (1997). Those variations lead to irregular cells and to non-periodic arrangement of the cell walls. Therefore, a quantitative study to investigate how the micro-structural variability can affect the macroscopic effective mechanical properties has been performed as a preliminary step. Statistical variations in the underlying models have been accounted for (Silva et al., 1995; Silva and Gibson, 1997; Zhu et al., 2001; Alkhader and Vural, 2008). Several methods described in Kraynik et al. (1991) account for a variability in the arrangement of cell walls of hexagonal honeycombs by modifying the initial two-dimensional unit cell analysis, Warren and Kraynik (1987). Those authors develop structure-property relationships for arrays of hexagonal cells endowed with varying sizes and shapes, but they conserve an angle of 120◦ between the three struts common to each node. The results of those authors lead to the conclusion that the specific spatial arrangement and size distribution of the unit cells hardly affect their elastic response. In order to generate a microstructural irregularity (or nonperiodicity), a spatial perturbation has been applied to the vertices of a regular triangular truss network in random directions (Der Burg et al., 1997; Chen et al., 1999; Chen and Fleck, 2002; Alkhader and Vural, 2008), expressed by the following equations: χ′i = χi + λr cos(θ)

γi′ = γi + λr sin(θ)

(1)

where θ is a uniformly distributed random variable, r is a random variable and χ′i and γi′ are the perturbed coordinates, with the non primed component being the original coordinates. λ is the perturbation parameter which specifies the degree of irregularity; note that it has been chosen in a manner that cell convexity is preserved. Fig. 1 shows different cellular structures generated for different values of λ. Each topology is 9 × 9

J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S

1835

4 (2011) 1833–1845

fact = 0.25

(a)

(b)

fact = 0.65

fact = 0.5

(c)

(d)

Fig. 1 – Cellular structures generated for FE analyses including stochastic triangular topologies. Sample microstructures subjected to various amounts of perturbation are shown. Note that nodal connectivity is 6 for these topologies. (a) Regular triangular topology. (b) The perturbation factor is 0.25, (c) 0.5, (d) 0.65.

cells and each beam element is 1 µm long within the initially regular topology. Finite element analyses were performed (using the FE software Abaqus implicit) in order to get the effective tensile moduli, which are then compared with their counterpart for the regular topology. Elastic linear and quasistatic frameworks have been considered; each cell element is modeled by three linear Bernoulli elements (Abaqus element type B23) considering a cubic formulation. Simulation results show that the perturbed topology in the structure introduced a small amount of variances in the elastic constants for an isotropic triangular unit cell (the coefficient of variance lies in the interval [0.2%–4%]. This finding is in agreement with many works in the literature, especially those dealing with cellular structural materials such as foam, Silva et al. (1995), Silva and Gibson (1997) and Zhu et al. (2001). The relative variation is represented versus the perturbation parameter λ on Fig. 2, showing that the non-periodic arrangement of cells does only weakly affect the elastic properties of the overall network, hence it is legitimate to adopt in the sequel a quasi periodicity assumption.

4 3.5 3 2.5 2 1.5 1 0.5 0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

Fig. 2 – Relative error in percent of the Young modulus versus the perturbation factor λ.

3.

Method and technique

One may classify the networks according to their connectivity: Three fold, four fold, and six fold connectivity

1836

J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S

networks are the most encountered types of biological membranes, Boal (2002), corresponding to the square, the hexagonal and the triangular networks respectively. The beams of the cell lattice represent the network chains, which are endowed with tensile and bending moduli evaluated from geometrical and mechanical parameters of the macromolecular chains to be given later.

3.1.

beam deformed

Discrete Homogenization technique

The discrete homogenization method is considered as a mathematical technique to derive the equivalent continuous medium behavior of a repetitive discrete structure made of elementary cells. This technique is inspired from the homogenization of periodic media developed thirty years ago by Sanchez, 1980, Bakhvalov and Panasenko, 1984, Panasenko 1983 and more recently applied by Warren and Byskov (2002) and Mourad (2003). It has been also combined with the energy method by Pradel and Sab and applied to discrete homogenization. More details about the method can be found in Caillerie et al. (2006); Mourad (2003). The developed method in Caillerie et al. (2006) and Mourad (2003) uses the moment equilibrium equation to replace the expression of the transverse forces, which is not straightforward within the considered application field, namely molecular dynamics or interatomic physics. Previous authors adopt the mechanics of interacting bars, whereas we choose beam mechanics for the description of the lattice behavior. The discrete homogenization method consists of assuming asymptotic series expansions of both the node displacements, tension and external forces in the successive powers of a small parameter labeled ε, defined as the ratio of a characteristic length of the basic cell to a characteristic length of the lattice structure. Those expansions are then inserted into the equilibrium equation, conveniently expressed in weak form. The balance equations of the nodes, forces–displacement relations and the moment–rotation relations of the beams are developed by inserting those series expansions and by using Taylor’s expansion of finite differences. The discrete sums are finally converted in the limit of a continuous density of beams into Riemann integrals, thereby highlighting continuous stress and strain measures. The calculations have been completed for a quite general truss and the results give a general and closed form expression of elastic properties. The method has given rise to implementation into a dedicated software, and can handle more complicated membranes. Note that second order development of the variables versus ε has been taken into account in order to capture the micropolar effect as mentioned in Warren and Byskov (2002), which is the main novelty of the method, especially for biological membranes.

3.2. Asymptotic parameters and description of the lattice geometry The discrete asymptotic technique requires the development of all variables as Taylor series vs. the small parameter ε, namely the beam length lεb , the beam width tεb , the lattice thickness eεb , the displacement uεn and the rotation at the lattices nodes φεn (they constitute the kinematic variables).

4 (2011) 1833–1845

med

efor und

beam

Fig. 3 – Kinematic and static parameters of a lattice beam.

The Bernoulli beam model is considered in this work. From the results of Mourad (2003), one can express the beam length as follows lb = lb0 + εlb1 + ε2 lb2 + · · · + εp lbp . The normal and transverse efforts, and the moment at the beam extremities can be successively expressed versus the kinematical nodal variables as, Dos Reis and Ganghoffer (2010): Nε =

Es S

(1Ubε · eb )   12Es Iz lb Ttε = (1Ub · eb⊥ ) − ε (φO(b)ε + φE(b)ε ) 2 (lb )3   12Es Iz 2 (lb )2 1 ε MO(b)ε = 2φO(b)ε + φE(b)ε − ε(∆bε · eb⊥ ) 6 2 (lb )3   b 2 12Es Iz 2 (l ) 1 ME(b)ε = φO(b)ε + 2φE(b)ε − ε(1Ubε · eb⊥ ) ε 6 2 (lb )3 lb

(2)

where Iz is the quadratic moment of the beam, eb the unit director for each beam and eb⊥ the transverse unit vector. MO(b)ε and ME(b)ε are the moment at the origin and the end positions of a generic beam respectively. The truss under consideration is made of beams and is completely defined by the positions of the nodes and their connectivity. Each beam links two nodes and is oriented so that it has an origin ˜ and an end node E(b). ˜ Although we can choose the node O(b) origin node as part of the reference cell, this is not necessarily the case for the end node, which nevertheless belongs to the next neighboring cell. Moreover, we associate to each beam a ˜ (Fig. 3). Each extremity node defined as its center denoted C(b) node has two displacements in the two principal directions and one rotation in the plane (i, j). The beam kinematic parameters together with the efforts and moments are shown in Fig. 3. We give the essentials relative to the asymptotic homogenization technique considering the micropolar framework, which is the more complete equivalent continuum, as it incorporates a microrotation in addition to the displacement as kinematic descriptors at the continuum level. The equilibrium of forces for the whole lattice writes in virtual power form and after asymptotic development as − − νi ∈Z2 b∈BR

Tεb (vε (O(b)) − vε (E(b))) = 0

(3)

4 (2011) 1833–1845

J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S

with v(.) a virtual velocity field choosing to vanish on the edges. The vector of effort Tb decomposes into a normal and a transverse contribution as Tb = Nb eb + Ttb eb⊥

(4)

with the normal and transverse efforts previously given in Eq. (2). The director vector eb (unit length, Fig. 3) and the beam lengths remain fixed under the adopted small strain framework. The asymptotic development of the virtual velocity and rotation rate are next expressed: For any virtual velocity field vε (λ), a Taylor series development leads to vε (O(b)) − vε (E(b)) = vε (λε + εδib ) − vε (λε ) ∂v(λε ) ib δ (5) ∂λi wherein we have parameterized any point within the surface element representative of the membrane by curvilinear coordinates λi (they may be conceived as Lagrangian coordinates of the material points); this allows treating membranes which exhibit a local curvature in their reference state. The discrete equilibrium of moments expresses in virtual power form as the following sum over all lattice nodes − − MO(b) w(O(b)) = ε

νi ∈Z2 b∈BR

+ME(b) w(E(b)) + lb (eb ∧ Ttb )w(C(b)) = 0

(6)

where w(C(b)) is the rotation rate of the incorporated central node of the beam, see Fig. 3. The rotation rate is similarly expanded taking into account the central node of the beam wO(b)ε = w(λ) 1 O(b)ε (w + wE(b)ε ) (7) 2 ∂w(λ) δi . wE(b)ε (λ + εδi ) = w(λ) + ε ∂λi In the forthcoming development, and for simplicity, a rectangular section of the beam will be considered, with a constant thickness e = 1 and a width tb . Thus, we can simplify the expression of the bending and stretching stiffnesses by defining

wC(b)ε =

b

the slenderness parameter η = tb (which has a finite value). l Inserting Eqs. (2) and (5) into (3) and considering the aforementioned simplifications, one obtains after some developments and ordering following the successive powers of ε the equilibrium of forces as Eq. (8)  − −[  ε2 Es η(1U1 · eb )eb + Es η3 1U1 · eb⊥ νi ∈Z2 b∈B

  1 O (b) E (b) Es η3 lb (φ0 R + φ0R ) eb⊥ 2   1 +ε3 Es η(1U2 · eb )eb + Es η3 1U2 · eb⊥ − Es η3 lb 2   ] ε) ∂φ ∂v(λ O (b) E (b) × φ1 R + φ1R + 0 δib eb⊥ · δi = 0 ∂λi ∂λi −

(8)

with 1U1 the first order difference of the displacement obtained versus ε as   ∂u E (b) O (b) 1Ubε = ε u1R − u1 R + 0 δib ∂λi    1U1

E (b) O (b) + ε2 (u2R − u2 R ) .





1U2



1837

The previous discrete equation can be transformed into a continuum Riemann integral on the (surface) domain Ω when the small parameter ε tends to zero: For any sufficient reg∑ ular function, the quantity ε2 νi ∈Z2 g(ενi ) can be interpreted  as the Riemann sum of an integral over Ω , Ω g(λ)dλ when ε → 0. Thus, the equilibrium equation in translation becomes after homogenization ∫ ∂v Si · dλ = 0 (10) ∂λi Ω evidencing a stress vector Si , which splits into a first and a second order contribution, viz. Si = Si1 + εSi2 , with Eq. (11).  − Si1 = Es η(1U1 · eb )eb + Es η3 1U1 · eb⊥ b∈BR

 1 O (b) E (b) Es η3 lb (φ0 R + φ0R ) eb⊥ 2  − Si2 = Es η(1U2 · eb )eb + Es η3 1U2 · eb⊥ b∈BR   1 ∂φ E (b) O (b) − Es η3 lb φ1 R + φ1R + 0 δib eb⊥ . 2 ∂λi −

(11)

Similarly to previous developments, the moment equilibrium (6) is homogenized, inserting the asymptotic expansion (7) of the virtual rotation rate. After simplifications and passing to the limit ε → 0 in the discrete sum, the moment equilibrium equation after homogenization takes the form ∫ ∂w µi · dλ = 0 (12) ∂λi Ω with the couple stress vector µi also identified on two orders, viz. µi = µi1 + εµi2 , with µi1 =

−

Es η3

(lb )3 ER (b) O (b) − φ0 R )δib (φ0 12

Es η3

(lb )3 12

b∈B

µi2 =

− b∈B



E (b)

φ1 R

O (b)

− φ1 R

+

 ∂φ0 ib ib δ δ . ∂λi

(13)

The very expression of the stress and couple stress vectors Si and µi in (10) and (12) is dependent on the unit cell topology and on the mechanical properties of the network beams. Observe that the homogenized equilibrium equations (10) and (12) involve virtual velocities (v, w), respectively equivalent to ˙ ˙ φ). the rates (u, In a second step, the equilibrium equations of the equivalent micropolar continuum are written in virtual power form, in order to highlight the stress tensor σ and the tensor of micromoment m as dyadic products of the evidenced stress and couple stress vectors with the gradient of the position vector with respect to the curvilinear coordinates. The following transformations from the Cartesian to the curvilinear coordinates λi are expressed, with R the position vector of any material point: ∂v ∂R = ▽x v · ; ∂λi ∂λi

∂w ∂R = ▽x w · ∂λi ∂λi

(14)

leading to the following equilibrium equation of the equivalent micropolar continuum ∫ ∫ ∂v ∂w (gσ · eiλ ) · dλ + (gm · eiλ ) · dλ = 0 Ω    ∂λi Ω    ∂λi Si

(9)

∫ → Ω

Si ·

µi

∂v ∂λi

∫ dλ + Ω

µi ·

∂w ∂λi

dλ = 0

(15)

1838

J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S

with g the Jacobian of the transformation from Cartesian to curvilinear coordinates. Hence, the last equation shows that the homogenized network obeys the same equilibrium equation as the equivalent continuum micropolar, with the stress and couple stress tensor identified with the following dyadic products: σ=

1 i ∂R ; S ⊗ g ∂λi

m=

1 i ∂R . µ ⊗ g ∂λi

(17) {m} = [C]{ε} + [D]{κ}. This form of the continuum constitutive law can presently be identified from the expressions of the homogenized stress and couple stress tensors together with Si and µi expressions: 1 1 ∂R ∂R + ε S2 ⊗ S ⊗ g 1 ∂λi g ∂λi       [A]{ε}

[B]{κ}

(18)

[D]{κ}

[C]{ε}

We restrict ourselves to centro-symmetric unit cells in this work, which entails the vanishing of the coupling matrices [B] and [C], as shown in Trovalusci and Masiani (1999). The expression of the stress and couple stress vectors substantially simplifies as follows, see Eq. (19) Si = Si1 =

Es η(1U1 · eb )eb +

−

−  Es η3 1U1 · eb⊥

b∈BR  1 OR (b) ER (b) 3 b − Es η l (φ0 + φ0 ) eb⊥ −2 b eb⊥ )δib = (Nb1 eb + Tt1

b∈BR

(19)

b∈B

µi = µi2 =

− b∈B

  ∂φ (lb )3 E (b) O (b) φ1R − φ1 R + 0 δib δib Es η3 12 ∂λi

− 1 E (b) O (b) = (M R − M2 R )δib 2 2 b∈BR

b , and Mn , respectively, the first order normal and with Nb1 , Tt1 2 transverse effort and the second order moment, obtained when expanding the expressions (2) versus the small parameter ε. Those expressions still involve the unknown displacements un1 , un2 and rotations φn0 , φn1 , which are determined for all nodes using the equilibrium equations (3) and (6).

4.

peptides

1.3 nm

sugar rings

b

{σ} = [A]{ε} + [B]{κ}

1 ∂R ∂R 1 m = µ1 ⊗ . + ε µ2 ⊗ g g ∂λi ∂λi      

2 nm

b

(16)

Using the symmetry properties of the lattice, we can simplify the expressions of the Cauchy stress σ and the couple stress m. The general form of the constitutive equations of linear micropolar elasticity relating the stress and couple stress tensors to the strain and curvature tensors is as follows

σ=

a

4 (2011) 1833–1845

Results and discussion

We shall consider two important types of biological membranes and evaluate their effective properties from the discrete homogenization technique exposed in the previous section.

θ

a

Fig. 4 – Peptidoglycan network a) dimensions of network chains (peptides and sugar rings) Koch and Woeste (1992), (b) face view of the section of the peptidoglycan network and definition of the geometrical model.

4.1. Equivalent mechanical properties of the peptidoglycan cell wall The peptidoglycan network is built from two non-equivalent chains, as shown in the face view of a section (Fig. 4) (from Koch and Woeste (1992)): Sugar rings run in one direction and peptide strings form transverse links. It is believed that this molecular anisotropic organization is dictated by design principle of the cell, such that the stiffer chains act as reinforcement in the direction that shall sustain the maximal stress. Note that the knowledge of the mechanical properties of peptidoglycan is of importance for understanding bacterial growth and form. The peptide and the glycan are modeled as beams with a regular circular cross section with radii respectively equal to 0.5 and 1 nm, Boal (2002); the angle θ is used as a descriptor of the topology of the glycan network (Fig. 5); it varies between 5◦ and 20◦ , according to the dimensions of the molecular chains of the unit cell shown in Fig. 4. The average Young moduli of glycan at low temperature (T = 273 K), is given by Ygly = 4.799 × 107 J/m3 and Ypep = 1.53 × 107 J/m3 for peptides: The Young moduli are elaborated from the persistence length according to the relation ξp = KEIzT (Boal, 2002), wherein E, Iz , T, KB are B the tensile modulus, the quadratic moment (dependent on the beam cross section), the absolute temperature and Boltzmann constant respectively. Considering the properties of individual chain, the peptide is endowed with the classical properties of entropic springs, given that its end-to-end length in the network lpep = ree = 1.3 nm is less than its contour length Lc = 4.2 nm. The persistence length of the peptide string can be obtained from the simplified relation ⟨r2ee ⟩ = 2ξp,peptide Lc , (the bracket denotes the average value) giving ξp,peptide = 0.2 nm. The glycan chains are comparatively much stiffer, with ξp,glycane ≃ 10 nm (Stokke and Brant, 1990) and lgly = ree = 2 nm. The whole lattice is generated from the repetition of the unit cell shown in Fig. 4(b) thanks to two periodicity vectors defined in the Cartesian basis. After calculations, we extract the homogenized moduli of the hexagonal lattice from the equivalent stiffness matrix, expressed versus

J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S

1839

4 (2011) 1833–1845

a

a 3 × 107

0.6

2.5 × 107

0.5

2 × 107

0.4

E1 E2

1.5 × 107

0.3

0.2

1 × 107

0.1

5 × 106

5

10

15

20

25

5

30

10

15

20

25

30

θ

θ

b

b 2.8 × 106 – 30 2.6

× 106

2.4

× 106

– 20

0

– 10

20

10

30

θ –1

G –2 2.2 × 106

2 × 106

–3

5

10

15

20

25

30

θ

–4

Fig. 5 – (a) Evolution of effective tensile moduli versus the geometrical parameter θ for the peptidoglycan. (b) Shear modulus.

Fig. 6 – Equivalent Poisson’s ratio versus the geometrical parameter θ for peptidoglycan network.

the geometrical and mechanical lattice parameters. Four mechanical parameters characterize the derived effective planar micropolar model: Two tensile moduli E1 ∗, E2 ∗, one shear modulus G∗, one Poisson’s ratio ν12 or ν21 , as the structure is shown to be orthotropic for all values of θ (the effective stiffness matrix is symmetrical), and two additional (micropolar) moduli γ∗ and κ∗. Those properties are next expressed versus the Young modulus, the slenderness ratio η and the configuration variable θ: [Eq. (20) is given in Box I.] We plot the evolution of the homogenized elastic properties versus the geometrical parameter θ in Figs. 6 and 7, with a range of values chosen in the interval [5◦ –30◦ ]. The equivalent tensile modulus along the x direction exhibits a maximum for a square topology (θ is nil) and shows strong variations around this maximum (it is an even function of θ). The shear modulus increases continuously in the selected range of variation of θ. The calculated effective properties appear to be very close to those found experimentally and to those calculated by other approaches. In Boal (2002), the shear modulus is about G∗ = 3.6 · 106 J/m2 , close to the homogenized value G∗ = 2 − 3 · 106 J/m2 . Moreover, the peptidoglycan has

two different effective Young moduli in the two principal directions, respectively E1 ∗ = 1 − 3.2 · 107 J/m3 and E2 ∗ = 3 − 5 · 106 J/m3 . The homogenized value is further close to measurements (Boal, 2002), which give E∗ = 2 − 3 · 107 J/m3 . A (symbolic) calculation based on the effective properties shows that the network is orthotropic (viewed as a continuum), which sets a relation between the two tensile moduli and the two Poisson’s ratios. The contraction coefficient ν21 increases monotonously versus θ, since the networks contracts more when the bars align in the y direction. The Poisson’s ratio ν12 shows a complex evolution pictured in Fig. 6, increasing through a maximum (the transverse contraction is maximal for an angle around 20◦ ) and decreasing thereafter. Considering a fixed geometry (the angle θ is fixed) and eliminating the variable slenderness ratio η, the scaling behavior of the equivalent elastic moduli versus the effective density is obtained, showing a nonlinear monotonous increase (Fig. 7). In order to assess the relevance of a micropolar continuum model, a micropolar characteristic length in bending lchar is evaluated versus the beam length, the angle θ and the slenderness ratio η in Eq. (21). This relation is extracted from

1840

J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S

E1 ∗ = E2 ∗ =

G∗ =

ν12 = ν21 =

K=

4 (2011) 1833–1845

η3 lgly Ygly cos(θ) (1 − cos2 (θ) + η2 cos2 (θ))(lpep + lgly sin(θ)) η3 Ygly Ypep (lpep + lgly sin(θ)) cos(θ)(2Ygly

η2 l

2 2 2 pep + lgly Ygly cos (θ) + η lgly Ypep (1 − cos (θ)))

η3 Ygly Ypep l2gly (lpep + lgly sin(θ)) cos(θ) (2η2 l2gly Ypep sin(θ)lpep + 2Ygly cos2 (θ)l3pep + l2pep Ypep lgly cos2 (θ) + η2 l3gly Ypep + l2pep η2 lgly Ypep − l2pep η2 lgly Ypep cos2 (θ)) lgly (η2 − 1) sin(θ) cos(θ) (1 − η2 cos2 (θ) − cos2 (θ))(lpep + lgly cos(θ)) (η2 − 1) sin(θ)(lpep + lgly sin(θ)) 2 2η Ygly lpep + lgly Ypep cos2 (θ) + lgly η2 Ypep cos2 (θ)

(lpep cos2 (θ) + lpep η2 sin2 (θ) + η2 lgly sin(θ) cos2 (θ) + η2 lgly sin3 (θ)η3 lgly 2 Ypep Ygly )

γ∗ =

(2 cos2 (θ)lpep 3 Ygly + η2 lgly 3 Ypep + 2Ygly η2 sin2 (θ)lpep 3 ) cos(θ) 3 3 1 cos(θ)Ygly η lgly 12 lpep + lgly sin(θ)

(20) Box I. the stiffness matrix (considering an orthotropic behavior)

a

l2char =

3 × 107

γ∗ 2(2µ ∗ +κ∗)

for a given direction (the moduli are extracted from the equivalent stiffness matrix), for the x-direction. 2 × 107

l2xchar =

1 × 107

0 0.1

0.3

0.5

0.7

0.9

ρ

b 4 × 106

1 b 2 η2 − 3 sin(θ) + η2 sin(θ) + 3 (l ) . 48 1 + sin(θ)

(21)

Two possible choices of the equivalent continuum have been considered, namely the Cauchy continuum and the micropolar continuum (the first one being a subcase of the latter), including additional rotational degrees of freedom compared to the former continuum model. The choice of the more appropriate continuum seems difficult a priori, since the equivalent properties may show unusual values as highlighted by previous results. Hence, we rather elaborate an a posteriori criterion for the choice of the equivalent continuum, relying on the comparison of the characteristic length of the micropolar continuum lchar to the beam length lb and to a macroscopic length at the scale of the lattice L:

3 × 106

- when lchar ≤ lb or lchar ≪ L, micropolar effects are unnoticeable at the macroscopic scale, and one may adopt a Cauchy continuum;

× 106

- when lb ≤ lchar ≤ L or when lchar ≈ L, micropolar effects have an impact at the macroscopic scale and it is justified to use a micropolar effective continuum.

2

1 × 106

0 0.1

0.3

0.5

0.7

0.9

ρ

Fig. 7 – Variation of the elastic moduli versus the density for different values of θ and different slenderness ratio η ∈ [0.01 − 0.5].

The ratio of the characteristic micropolar length lxchar (Eq. (21)) to the characteristic unit cell size 2lb cos(θ) is plotted for the peptidoglycan network, with the previously adopted mechanical parameters (Fig. 8). For classical configurations of the network (with positive θ), the ratio is less than unity, showing that micropolar effects may be neglected (they act at a length scale lower than the unit cell size). For strongly re-entrant hexagonal lattices corresponding to angles θ ≤ −35◦ , the ratio is larger than unity, thereby meaning that the microrotation effect becomes noticeable at a macroscopic scale, hence a Cosserat medium

J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S

1841

4 (2011) 1833–1845

a

14

12

10

8

6

b

4

2

θ – 40

–30

–20

–10

0

10

20

θ

Fig. 8 – Evolution of the ratio of the characteristic micropolar length to the unit cell size of the peptidoglycan network versus θ. is appropriate for such networks. A negative θ (corresponding to a so-called re-entrant lattice) entails a negative Poisson’s ratio (Fig. 6(b)), although it is not related to micropolar effects, with values becoming as negative as −5 when θ reaches − π3 . The mechanism responsible for the negative Poisson’s ratio is the deployment of the re-entrant beams in the direction transverse to the applied loading. In addition to the persistence length used to describe the correlation of chain segments orientations at the molecular scale, an additional length at the microscopic scale is introduced, the micropolar length, describing flexion effects at the scale of the unit cell. The present study shows that micropolar effects naturally emerge as an outcome of the homogenization procedure, and the micropolar length may have an impact at the scale of the whole membrane for certain topologies of the filament network. Biological materials such as bone were already reported to exhibit size effects amenable to a Cosserat model (Lakes, 1995), slender specimens having a higher apparent stiffness than thick ones. We presently show that biological networks may also exhibit such size effects and provide a methodology based on multiscale modeling to quantify those effects. The rationale for introducing micropolar effects in biological networks is that the filaments in the network are generally convoluted (the end to end distance may be much smaller than the contour length); hence, bending and twisting due to local rotations may be the dominant deformation modes of those networks. The idea of a couple stress can be traced to Voigt during the early development of the theory of elasticity; the concept of micropolar solid was fully developed one century ago by the Cosserat brothers (1909), but it remained a purely conceptual object until man was able to produce and characterize solids (either artificial of natural) exhibiting such scale effects. Later on, theories incorporating couple stresses were developed using the full capabilities of modern continuum mechanics, (Mindlin, 1965). As a consequence of the micropolar behavior, a size effect appears in torsion and bending, but not in tension; the Poisson’s ratio is unaffected.

Fig. 9 – Triangular network of the erythrocyte cytoskeleton (a) with its idealization configuration parametrized by the angle θ (b).

4.2. Equivalent mechanical properties of the Erythrocytes network The discrete homogenization technique is further exemplified for the erythrocyte (red blood cell, abbreviated RBC in the sequel) membrane. The diversity of models for the mechanics of RBCs from the literature can be divided into microscopic molecular based models (discrete models) and macroscopic continuum models, Hartmann (2010) and references therein. Among them, one can distinguish static models based on energy minimization over the whole cell (Discher et al., 1998), or dynamic models (Noguchi and Gompper, 2005), including the fluid flow around the RBC. The membrane of the erythrocyte has both a low density and a high flexibility. Furthermore, the corresponding lattice has six-fold connectivity, with however 3% and 8% of the junctions presenting five or seven-folds respectively (Liu et al., 1987; Mohandas and Evan, 1994; Trovalusci and Masiani, 1999; Feyel and Chaboche, 2000). The description of the in-plane network, parameterized with θ, is given in Fig. 9. Considering a low temperature (close to T = 273 K) and an initially non-stressed membrane, the spring constant for any chain in the network is considered as constant, with an average value given in Boey et al. (1998) βkef f a2 = 0.23

(22)

with β the inverse temperature K1T , a = 6.4 nm the bead B diameter Boey et al. (1998) and kef f = YS the tensile modulus lb of a beam. Using Eq. (22), the Young modulus y for any chain in the network is then calculated from the following expression Y=

0.23 lb βa2 S

(23)

with S the section of the considered chain. In the sequel, we consider a rectangular beam section with a chain width t = 9.96·10−3 nm (hence assigning the quadratic moment) and an average end to-end beam length lb = 200 nm (Feyel and Chaboche, 2000). Each tetramer is characterized by a constant beam length and an angle θ varying in the neighborhood

1842

J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S

of the reference value θ = 60◦ , which allows deriving the equivalent properties for different configurations (Fig. 10). The normalized characteristic length of the triangular unit cell is in fact less than 3 · 10−6 for θ larger than 10◦ , hence a Cauchy continuum is sufficient for the triangular erythrocyte network, as the micropolar effects appear at a length scale much lower than the unit cell size. The equivalent properties for the selected homogenized Cauchy continuum are evaluated as Eq. (24) E1 ∗ = − E2 ∗ = − G∗ = ν12 = ν21 =

0.00005

0.00004

E1 E2

0.00003

ηY sin(θ)(η2 cos(θ) − cos(θ) + 1) (cos(θ) − 1)(η2 cos2 (θ) − cos2 (θ) + 1)

0.00002

ηY sin(θ)(η2 cos(θ) − cos(θ) + 1) cos(θ)(η2 cos2 (θ) − cos2 (θ) + cos(θ) − η2 cos(θ) − 1)

ηY sin(θ)(η2 − η2 cos(θ) + 4 cos3 (θ) + cos2 (θ)) sin(θ)(1 + 4 cos(θ))

50

(1 − η2 ) cos2 (θ)

(1 − η2 )(cos(θ) − cos2 (θ)) − 1

60

65

70

b

(1 − η2 ) cos2 (θ) + 1 (1 − η2 )(cos(θ) − 1) cos(θ)

55

θ

0.000011

.

(24)

The evolution of the effective tensile moduli versus θ is pictured in Fig. 10(a). In the particular case of an equilateral triangle (θ = 60◦ ), the elastic moduli and the Poisson’s ratio in the two principal directions are identical (both properties are equal) due to the fact the equivalent continuum material is then isotropic. The elastic modulus in the transverse direction E2∗ becomes weak when the triangle is flattened in the vertical direction (the angle θ is decreased), whereas the modulus E1∗ increases, as all threads tend to align in the direction. The Poisson’s ratio of the RBC membrane is close to 0.33 for the reference angle θ = 60◦ . Note that due to the small perturbation framework, consistent with low temperatures, the topological fluctuations of this network are not important enough to generate a negative Poisson ratio due to entropic effects, as reported in Boal et al. (1993). We can nevertheless notice that Poisson’s ratio tends to zero for angular deviations from the reference value of (60◦ ) about 10%. The shear modulus is plotted versus θ on Fig. 10(b) and decreases when θ increases (the triangle elongates and becomes flat in the y-direction). For an equilateral tetramer, the effective shear modulus is close to G∗ = 9.52 × 10−6 J/m2 , which agrees with experimental findings (Lenormand et al., 2001; Hartmann, 2010; Mohandas and Evan, 1994). Dao et al. (2006) find a uniaxial tension Young’s modulus E ≈ 2.21 · 10−5 J/m3 , close to the calculated effective Young modulus E1∗ = E2∗ = 2.5 · 10−5 J/m3 . Moreover, the planar modulus of compressibility takes the value KA = 14.6 · 10−6 J/m2 and agrees well with values from the literature (Boal, 2002; Hartmann, 2010).

4.3.

a

4 (2011) 1833–1845

Insight into the nonlinear behavior of membranes

As biological membranes undergo thermal fluctuations, they are prone to undergoing large configurational changes: In order to capture the nonlinear mechanical behavior of biological membranes on the basis of the calculated effective linear elastic behavior (the lattice geometry is fixed), a perturbative method is adopted, aiming at a first modeling approach of the geometrical nonlinearity due to the change

0.000010

0.000009

G

0.000008

45

50

55

60

65

70

θ

Fig. 10 – (a) Evolution of the tensile moduli versus θ. (b) Shear modulus.

of lattice configuration. We presently do not model the chain of events leading to such a modification, but limit ourselves to the mechanical factors responsible for such configurational changes. The biochemical mechanisms responsible for shape changes of the membrane have been evoked in Section 2, in relation to remodeling. The present model does not describe remodeling, a process relaxing the free energy, but focuses on deformation induced by pure mechanical forces. The nonlinear elastic behavior of some biological membranes including the RBC has been acknowledged in several contributions; instead of adjusting in an ad-hoc manner an hyperelastic potential such as in Gong et al. (2009), we shall build the nonlinear response of membranes in a numerical manner from discrete homogenization. We shall accordingly extend the linear framework so far adopted and consider the impact of a variation of the lattice geometry on the equivalent moduli, and hence on the effective membrane behavior. The adopted perturbative scheme relies on the following steps: i. Fix a reference network defined by one or a few geometrical descriptors; the equivalent moduli are known for this chosen reference configuration.

J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S

3

1843

4 (2011) 1833–1845

× 106

σ 11 (Pa)

2.5

2

1.5

1

0.5 0.02

ε

ε 11 0.04

0.06

0.08

0.1

0.12

0.14

0.16

Fig. 11 – Stress strain response of the peptidoglycan network.

Fig. 12 – Evolution versus strain of the planar compressibility modulus of the erythrocyte.

ii. Perturb the network geometry by allowing a small variation of the geometrical parameters chosen as descriptors, in a displacement control approach. Evaluate then the new updated equivalent moduli for this new geometry. iii. Elaborate strain measures from the geometrical parameters and record the variation of the moduli versus those strain measures. The functional dependence of the updated moduli is deemed representative of the nonlinear mechanical behavior of the membrane.

4.4.

The previous cases of the hexagonal and triangular lattices are successively considered, with a specific exploitation in terms of nonlinear properties for each lattice. The equivalent tensile moduli of the hexagonal lattices are calculated versus the strain variable for the hexagonal lattice, defined as the relative length ε=

l − l0 l0

(25)

with the actual length given by l = 2lb cos(θ), and an initial length l0 = 2lb cos(θ0 ), obtained from the previous expression for the reference angle θ = π/6. Accordingly, the tensile behavior of this network in the x-direction is simulated (at the nano-level) by a continuous decrease of the configuration angle θ, inducing in turn a continuous variation of the associated strain following (25). This is realized by controlling the displacement of the vertical right edge of the unit cell. The nonlinear response shows a geometrical hardening (the slope is increasing) in the xdirection traducing the progressive orientation of the threads into the x-direction when the unit cell length is increased (see Fig. 11). A kinematic locking will arise for a certain strain value (about 15% in strain). In the same vein, the planar compressibility modulus tr(σ) KA = δA A = tr(ε) , ratio of the trace of the stress to the trace of the strain (traducing the change of the local area), is evaluated for the triangular network representative of the erythrocyte, with geometrical nonlinearities captured by the same iterative procedure. It decreases by a factor two over a range of strains up to 10% Fig. 12, showing a tendency of the erythrocyte network area to become more expandable for large strains; this partly explains the ease of deformation of erythrocytes at large strains (although shear may be the dominant deformation mode).

Discussion and limitations

Although our results hold promise to reliably predict the mechanical behavior of biomembranes up to the nonlinear regime, the present modeling approach suffers from some limitations which we next underline. The membrane has here been treated as a single layer, which is a strong approximation as to the in-plane shear behavior (the in-plane stretch is much less influenced), since biological membranes consist of two leaflets endowed with a relative sliding mobility. This is however a common viewpoint in the literature, since the molecular based models treat the spectrin network and lipid membrane (which rather has a fluid like behavior) in a globalized manner (Dao et al., 2006). The interest of the developed homogenization technique lies in its flexibility and capability to handle any planar lattice, with the molecular chains assimilated to beams. In the elastic case, and for small perturbations, the effective properties are derivable as closed-form expressions of the geometry of the molecular chains (characterized by the persistence length, the quadratic moment) and mechanical properties (Young modulus). The persistence length is the geometrical parameter at the nanoscale that may be used as an input to account for the entropic elasticity of the network chains. The beam length in the network model may be assimilated to the average end-to-end distance of a given chain, which is known to fluctuate with temperature and depends upon the persistence length (Boal, 2002). This effect is only indirectly incorporated into the model, by changing the beam length according to temperature and the estimated value of the persistence length. The role of fluctuations of the chain segments are however not considered in the present model. The role of the fluctuations of the geometry and properties (the Young modulus is linked to the persistence length, itself affected by the temperature) over long distances may be assessed from a sensitivity analysis using the closed form expressions of the effective moduli. Developing this idea, one may even evaluate the nanoscopic parameters of a given membrane in a kind of inverse approach from measurements performed at a more macroscopic scale, using the expressions of the moduli to make such an adjustment. The network’s chains have been modeled as linear elastic (characterized by a tensile and a bending modulus), which may be considered as the small

1844

J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S

strain approximation of nonlinear models, such as the wormlike chain model. The further interest of the homogenization method is that it delivers the full compliance (or rigidity) matrix, reflecting the sometimes complex anisotropy of the sobuilt equivalent continuum (the material symmetry group of the discrete lattice is included in the material symmetry group of the homogenized membrane Trovalusci and Masiani (1999)). The two considered effective media, of Cauchy or Cosserat type, provide the same tensile moduli. The Cosserat media provides further micropolar moduli, and an associated characteristic length, the importance of which deciding upon the scale at which micropolar effects are noticeable. The studied membranes lead in their regular configurations (for low temperatures) to a classical Cauchy continuum; as an original aspect of this work, the reentrant hexagonal network of the peptidoglycan shows a characteristic length larger than the size of the unit cell, hence micropolar effects show an impact at the scale of the whole network, and have to be accounted for in the homogenized constitutive law. Due to micropolar effects, the shear modulus decreases µ∗ = µ − κ/2 and even becomes negative. Hence, biological networks will deform more easily in shear and flexion/torsion, especially in the large deformation regime. The Cosserat (micropolar) continuum model is therefore relevant for modeling membranes in specific configurations for which such exotic effects (negative contraction) become pronounced. Although the likelihood of such configurations may be low according to the induced energy cost and the resulting low probability (expressed by the factor exp(−1E/KB T)), they may locally be present due to thermal fluctuations. The extension of the present work to homogenization of membranes in both the geometrical and mechanical nonlinear frameworks (it is clear that topological fluctuations play an important role in the nonlinear behavior of membranes, in addition to the nonlinear mechanics of the chains) shall be addressed in the future, allowing the proper construction of hyperelastic potentials reflecting the membrane anisotropy in their very expression. Since the method is not based on an energy minimization principle but on the more general principle of virtual work, it allows a possible extension to viscous effects of the chains, an issue to be addressed in future contributions. Although planar lattices have been considered, the present homogenization technique allows modeling lattice shells (membranes endowed with a curvature in their natural state), provided suitable curvilinear coordinates are adopted. Contrary to models based on energy minimization, and as an intermediate mesoscopic scale is considered, the behavior is local and does not reflect the membrane curvature, which shall be taken into account in macroscopic models of the whole cell. Micromechanical analyses of biological networks based on those complex constitutive laws pave the way to multiscale finite element simulations over whole cells submitted to loadings reflecting the effect of their biological environment, or in order to analyse the deformation of the cell under a controlled loading. In particular, one could combine within a two scales calculation a homogenization–projection method to locally update the constitutive law together with the structural calculation over the whole membrane as in Feyel and Chaboche (2000).

5.

4 (2011) 1833–1845

Conclusion

A systematic micromechanical method for the calculation of the equivalent mechanical properties of biological networks has been developed, considering the membrane at the discrete level as the repetition of a unit cell, following the discrete asymptotic homogenization method. The beams of the cell lattice represent the threads, which are endowed with tensile and bending moduli evaluated from geometrical and mechanical parameters of the underlying macromolecular chains. Its knowledge together with Boltzmann factor and the quadratic moment of the beam allow evaluating the tensile and the bending moduli of each thread. Preliminary simulations performed over a perturbed lattice have shown that the assumption of a quasi periodic lattice may be reasonably adopted. The derived constitutive laws exhibit the anisotropy of the equivalent continuum membrane at a mesoscopic scale, which may be difficult to access experimentally or by other analytic techniques. The homogenized constitutive behavior delivers in the small perturbation case closed form expressions of the effective mechanical moduli versus the lattice topology and the mechanical properties of the thread, accounting for entropic elasticity of the threads and for an internal bending length responsible for micropolar effects, which may be significant at the continuum level. The nonlinear mechanical response of the membrane network has been obtained from a perturbation of the lattice geometry. REFERENCES

Alkhader, M., Vural, M., 2008. Mechanical response of cellular solids: role of cellular topology and microstructural irregularity. International Journal of Engineering Science 46, 1035–1051. Boal, D., 2002. Mechanics of the Cell. The Press Syndicate of the University of Cambridge. Boal, D., Seifert, U., Shillcock, J., 1993. Negative poisson ratio for two-dimensional networks under tension. Physical Review E 48 (6), 4274–4283. Boey, S., Boal, D., Discher, D., 1998. Simulations of the erythrocyte cytoskeleton at large deformation. Biophysical Journal 75, 1573–1583. Caillerie, D., Mourad, A., Raoult, A., 2006. Discrete homogenization in graphene sheet modeling. Journal of Elasticity 84, 33–68. Chen, C., Fleck, N., 2002. Size effects in the constrained deformation of metallic foam. Journal of the Mechanics and Physics of Solids 50, 955–977. Chen, C., Lu, T., Fleck, N., 1999. Effect of imperfections of the yielding of two dimensional foams. Journal of the Mechanics and Physics of Solids 47, 2235–2272. Dao, M., Li, J., Suresh, S., 2006. Molecularly based analysis of deformation of spectrin network and human erythrocyte. Material Sciences and Engineering C 26, 1232–1244. Der Burg, M., Shulmeister, V., Geissen, E., Marissen, R., 1997. On the linear elastic properties of regular and random open-cell foam models. Journal of Cellular Plastics 33, 31–54. Discher, D., Boal, D., Boey, S., 1998. Simulations of the erythrocyte cytoskeleton at large deformation II. Biophysical Journal 75, 1584–1597. Dos Reis, F, Ganghoffer, J -F, 2010. Discrete homogenization of architectured materials: Implementation of the method in a simulation tool for the systematic prediction of their effective elastic properties. Technische Mechanik 30, 85–109.

J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S

Feyel, F., Chaboche, J.-L., 2000. FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials. Computer Methods in Applied Mechanics and Engineering 183, 309–330. Gong, X., Sugiyama, K., Takagi, S., Matsumoto, Y., 2009. The deformation behavior of multiple red blood celss in a capillary vessel. Journal Biomechanical Engineering 131. Hansen, J., Skalak, R., Chien, S., Hoger, A., 1997. Influence of network topology on the elasticity of the red blood cell membrane skeleton. Biophysical Journal 72, 2369–2381. Hartmann, D., 2010. A multiscale model for red blood cell mechanics. Biomechanics and Modelling in Mechanobiology 9 (1), 1–17. Koch, A., Woeste, S., 1992. Elasticity of the sacculus of escherichia coli. Journal of Bacteriology 174, 327–341. Kozlov, M., McMahon, H., Chernomorik, L., 2010. Protein-driven membrane stresses in fusion and fission. Trend in Biochemical Sciences 35, 699–706. Kraynik, A., Reinelt, D., Princen, H., 1991. The nonlinear elastic behaviour of polydisperse hexagonal foams and concentrated emulsion. Journal of Rheology 35, 1235–1253. Lakes, R., 1995. On the torsional properties of single osteons. Journal of Biomechanics 28, 1409–1410. Lenormand, G., Hénon, S., Richert, A., Simeon, J., Gallet, F., 2001. Direct measurement of the area expansion and shear moduli of the human red blood cell membrane skeleton. Biophysical Journal 81, 43–56. Lim, C., Zhou, E., Quek, S., 2006. Mechanical models for living cells, a review. Journal of Biomechanics 39, 195–216. Liu, S.-C., Derick, L., Palek, J., 1987. Visualization of the hexagonal lattice in the erythrocyte membrane skeleton. Journal of cell Biology 104, 527–536. Mindlin, R., 1965. Stress functions for Cosserat continuum. International Journal of Solids and Structures 1, 265–271.

4 (2011) 1833–1845

1845

Mohandas, N., Evan, E., 1994. Mechanical properties of the red cell membrane in relation to molecular structure and genetic defects. Annual Review of biophysics and biomolecular structure 23, 787–818. Mourad, A., 2003. Description topologique de l’architecture fibreuse et modélisation mécanique du myocarde. Ph.D. Thesis, Institut National Polytechnique de Grenoble. Noguchi, H., Gompper, G., 2005. Shape transitions of fluid vesicles and red blood cells in capillary flows. Proceedings of the National Academy of Sciences of the United States of Amercia 102 (14), 14159–14164. Silva, M., Gibson, L., 1997. The effects of non-periodic microstructure and defects on the compressive strength of twodimensional cellular solid. International Journal of Mechanical Sciences 39, 549–563. Silva, M., Wilson, C., Gibson, L., 1995. The effect of non-periodic microstructure on the elastic properties of two dimensional cellular solids. International Journal of Mechanical Sciences 37, 1161–1177. Stokke, T., Brant, D., 1990. The reliability of worm-like polysaccharide chain dimension estimated from electron micrographs. Biopolymers 30, 1161–1181. Trovalusci, P., Masiani, R., 1999. Material symmetries of micropolar continua equivalent to lattices. International Journal of Solids and Structures 36, 2091–2108. Warren, W., Byskov, E., 2002. Three-fold symmetry restrictions on two-dimensional micropolar material. European Journal of Mechanics A/Solid 21, 779–792. Warren, W., Kraynik, A., 1987. Foam mechanics: the linear elastic response of two-dimensional spatially periodic cellular material. Mechanics of Materials 6, 27–37. Zhu, H., Hobdell, J., Windle, A., 2001. Effects of cell irregularity on the elastic properties of 2d Voronoi honeycombs. Journal of the Mechanics and Physics of Solids 49, 857–870.