## Abstract

Biological tissues possess the ability to adapt according to the respective local loading conditions, which results in growth and remodelling phenomena. The main goal of this work is the development of a new remodelling approach that, on the one hand, reflects the alignment of fibrous soft biological tissue with respect to representative loading directions. On the other hand, the continuum approach proposed is based on a sound micro-mechanically motivated formulation. To be specific, use of a worm-like chain model is made to describe the behaviour of long-chain molecules as present in, for instance, collageneous tissues. The extension of such a one-dimensional constitutive equation to the three-dimensional macroscopic level is performed by means of a microsphere formulation. Inherent with the algorithmic treatment of this type of modelling approach, a finite number of unit vectors is considered for the numerical integration over the domain of the unit sphere. As a key aspect of this contribution, remodelling is incorporated by setting up evolution equations for the referential orientations of these integration directions. Accordingly, the unit vectors considered now allow interpretation as internal variables, which characterize the material’s anisotropic properties. Several numerical studies underline the applicability of the model that, moreover, nicely fits into iterative finite element formulations so that general boundary value problems can be solved.

## 1. Introduction

Adaptation is a key biomechanical phenomenon occurring in hard as well as soft biological tissues. The modelling and simulation of such growth and remodelling processes constitute an active and continuously increasing field of research; for an overview, the reader is referred to the monographs by Fung (1990) and Humphrey (2002) as well as to the contributions in Holzapfel & Ogden (2003), the survey by Taber (1995) and the contribution by Cowin (2004). While various biological and chemical effects determine the macroscopic behaviour of the particular biological tissue of interest, we here restrict ourselves to a purely mechanical and non-coupled

approach. In general, biological tissues exhibit changes in mass and internal structure. This contribution exclusively focuses on the latter, which we denote as remodelling.

Fibrous soft biological tissues—for example, ligaments, tendons, muscles and skin—possess a pronounced composite-type multi-scale structure together with strongly direction-dependent or rather anisotropic mechanical properties (e.g. Silver *et al*. 2003). The local mechanical response of these material structures is typically determined by elastin and collagen fibres, respectively cable-like bundles thereof. Placing special emphasis on the latter, collagen molecules are secreted by the so-called fibroblasts, which order themselves as polymeric collagen fibrils. Based on these, cable-like collagen fibres are formed that turn out to be aligned with respect to the fibroblasts; see Alberts *et al*. (1994) for detailed background information. Influenced by genetic information as well as by local mechanical loading conditions, remodelling and self-assembly of these oriented fibre networks occur—the related processes often being denoted as fibre reorientation or rather turnover. Several experimental investigations on these remodelling effects are reported in the literature, even though the underlying mechanisms, based on which the respective cells sense the various mechanical and biochemical stimuli, until now remain fairly unknown. An overview on the adaptation of different types of biological tissue is provided in the survey by Wang & Thampatty (2006). To give two particular examples, adaptation of tendon-related tissue is discussed by Calve *et al*. (2004), whereas Breen (2000) studies tissue remodelling of the lung. The computational remodelling framework subsequently developed in this contribution is partly motivated by the investigations reported in Eastwood *et al*. (1998), in which a fibroblast-populated collagen lattice was tested. As a result, the macroscopically tension-type mechanical loads, including cell contraction, rendered the initially unstructured collagen fibre network to reorient with the local dominant stretch direction. In this regard, an anisotropic micromechanically motivated model that additionally incorporates this time-dependent remodelling effect is developed as this work proceeds.

Several continuum modelling approaches to describe the anisotropic material properties of soft biological tissues have been suggested in the literature; see, for instance, Almeida & Spilker (1998), the contributions in Cowin & Humphrey (2000), Alastrué *et al*. (2007) and Menzel & Steinmann (2001) for a comparison of two different modelling approaches. While phenomenological formulations are used to consider individual fibre families to be perfectly aligned with pre-defined directions, recent investigations additionally account for the dispersion of fibres; see Gasser *et al*. (2006) for a formulation based on generalized structural tensors and Marquez (2006) in view of experimental investigations to identify the orientation of fibre distributions. In this contribution, we will incorporate the anisotropic fibre contributions by means of a microsphere model. The computational modelling based on such a multi-scale-type unit-sphere ansatz dates back to the development of the so-called microplane models as embedded into small strain kinematics; for an overview, the reader is referred to Carol *et al*. (2001), Kuhl *et al*. (2001), Kuhl & Ramm (2000) and references cited in these works. The extension of the microplane framework to the finite deformation microsphere model was first developed for rubber elasticity, as well as rubber inelasticity, in the work documented in the series of papers by Miehe *et al*. (2004), Miehe & Göktepe (2005) and Göktepe & Miehe (2005). Recently, the model has also been applied to the simulation of elastic but anisotropic and pre-stressed blood vessels; see Alastrué *et al.* (2009*a*), in which the integration over the underlying unit sphere was additionally weighted by an orientation distribution function (odf) ρ with
To give a brief classification of some generic homogenization schemes, consider 〈•〉 to represent an appropriate averaging, respectively homogenization operator, that may include directional dependencies by means of ρ. Moreover, let ** S** denote a suitable stress measure, the computation of which may be based on
In order to capture the material’s anisotropy, a generalized structural tensor 〈

**⊗**

*e***〉 can be introduced, which together with the stretch quantity**

*e***determines the stresses. A remodelling formulation of approach (i) is proposed in Menzel**

*C**et al*. (2008). In view of further background information on the simulation of deformation-induced anisotropy evolution based on evolving structural tensors, the reader is referred to Menzel & Steinmann (2003

*a*,

*b*) and references cited therein. In this contribution, we make use of approach (ii), namely the microsphere model, which is extended to the incorporation of remodelling effects. The underlying anisotropy evolution is thereby directly linked to the numerical integration scheme applied to the unit sphere. At this stage, the well-established algorithm based on 42 integrations directions is applied (see Bažant & Oh 1986). Further integration schemes are provided in, for instance, Heo & Xu (2001), while their respective properties with application to anisotropic elasticity are discussed in Alastrué

*et al*. (2009

*a*). Additional modifications of such integration algorithms to efficiently recapture the material’s anisotropy are elaborated in Alastrué

*et al.*(2009

*b*) with emphasis on uniaxial odfs. Even though not elaborated in this work, both approaches (i) and (ii) can, in general, be extended to the incorporation of higher order moments such as 〈

**⊗**

*e***⊗**

*e***⊗**

*e***〉 and so forth. Furthermore, the evolution of the odf—in contrast to the direct local form of evolution equations used here—can alternatively be captured by means of global balance equations; see, for example, Papenfuss (2000) or the monograph by Virga (1994) with application to liquid crystals. Finally, a so-called FE**

*e*^{2}strategy could be adopted, wherein the macroscopic constitutive behaviour is determined from the solution of a sub-boundary value problem; see Miehe & Koch (2002) or Kouznetsova

*et al.*(2004). Such an approach (iii), however, is not in the focus of this work.

The continuum modelling of the adaptation of biological tissues is commonly either embedded into the framework of open-system thermodynamics or based on the theory of mixtures. A survey on open systems is provided in the monograph by Katchalsky & Curran (1965), while related continuum approaches with particular emphasis on growth phenomena have been addressed by, for example, Epstein & Maugin (2000), DiCarlo & Quiligotti (2002), Lubarda & Hoger (2002), Kuhl and Steinmann (2003*a*,*b*), Garikipati *et al.* (2004) and Guillou & Ogden (2006). In view of mixture-theory-based growth formulations, we refer the reader to the references by Humphrey & Rajagopal (2002), Quiligotti (2002) and Klisch & Hoger (2003). In this contribution, focus is placed on the remodelling of a single solid phase. The particular approach developed extends reorientation formulations previously proposed for the alignment of individual fibre directions as discussed in Imatani and Maugin (2002), Driessen *et al.* (2004) and Menzel (2005, 2007). In particular, we will make use of an evolution equation in terms of a shear-related driving force that reorients the respective direction vectors. Similar approaches have also been applied to the modelling of glassy polymers; see, for instance, Warren (1995) or Chen *et al.* (1996), in which macroscopic evolution equations were derived from odfs.

Long-chain molecules within fibrous soft biological tissues can, similar to rubber, appropriately be modelled within the framework of statistical mechanics; see the monograph by Flory (1969) for detailed background information. Here, we will refer the constitutive behaviour of a single chain to the so-called worm-like chain model; see Beatty (2004) for a review. While the isotropic chain network response can be determined by means of representative unit cells aligned with the principal stretch directions—a typical example being the eight-chain model advocated by Arruda & Boyce (1993)—anisotropic extensions reflecting transversely isotropic or orthotropic response have been suggested in the literature as well; see Bischoff *et al.* (2002) and Kuhl *et al.* (2006). As this work proceeds, the network will be captured by the microsphere framework, whereby we—apart from isotropic contributions—restrict ourselves to an affine fibre model. The extension to a non-affine formulation has been developed in Miehe *et al.* (2004) (see also the discussion by Chandran & Barocas (2006)), and could also can be adopted for the approach proposed in this work.

The paper is organized as follows: § 2 briefly reviews essential kinematic relations, based on which key aspects of the microsphere model are addressed in § 3. Section 4 constitutes the main body of this contribution, wherein the remodelling formulation is developed. Several numerical examples for homogeneous deformations and a non-homogeneous deformation are discussed in § 5, before the paper closes with a short summary in § 6.

## 2. Essential kinematics

Let describe the motion of the body considered, which transforms referential position vectors of material particles to their spatial counterpart . Moreover, the motion gradient and the right Cauchy–Green deformation tensor as well as their isochoric parts are denoted by standard notations, namely
2.1
with *J*=det(** F**)>0. In view of the subsequently elaborated computational homogenization scheme, we next introduce the kinematics referred to the underlying unit sphere . In this regard, let denote a unit vector or rather an element of the unit sphere, so that the affine stretches in the direction of

**are determined by the macroscopic deformation tensors via 2.2 Non-affine models, as based on additional representative deformation measures, have been proposed in the literature. At this stage, however, we restrict the unit sphere-based contributions to affine deformations.**

*e*## 3. Hyper-elastic microsphere model

Apart from the particular remodelling framework discussed later on, we make use of a hyper-elastic form so that the stress tensors of interest can be derived from a potential. In the following, we do not place emphasis on a full thermodynamical framework, which here should be embedded into the theory of open systems and additionally account for a so-called chemical potential, but solely outline the solid mechanical part of the continuum. Its contribution to the dissipation inequality results for the isothermal case in
3.1
wherein represents the material time derivative and Ψ denotes the Helmholtz free energy. As the material properties of the types of soft tissue we are interested in may be highly anisotropic, this energy function is assumed not only to depend on the deformation gradient ** F** but also on a finite number of direction vectors, namely with , so that Ψ(

**,**

*C*

*r*_{i};

**). Accordingly, the hyper-elastic form for the Piola–Kirchhoff stresses**

*X***renders the reduced dissipation inequality as 3.2 Apparently, corresponds to deformation-induced anisotropy evolution. As various soft biological tissues have been reported to respond quasi-incompressible, we adopt the well-established volumetric–isochoric split and make use of decomposing the isochoric part into an isotropic and an anisotropic contribution; to be specific 3.3 wherein possible heterogeneities of the material are not explicitly indicated.**

*S*The anisotropic part of the Helmholtz free energy, which allows incorporation of the fibrous properties of the biological tissue, is now further specified. In particular, a so-called microsphere formulation is adopted, so that the macroscopic stress quantity corresponding to Ψ^{ani} is determined by integration over the unit sphere. From the micromechanical point of view, the density of fibres, with respect to the direction, contributes to these macroscopic stresses. When assuming this collection of *n* fibres to be loaded or rather stretched in terms of the macroscopic isochoric deformation measure , i.e. we deal with an affine deformation, then Ψ^{ani} can be computed by means of the fibre-related Helmholtz free energy ψ^{ani} via
3.4
wherein is the respective direction considered for the integration over the unit sphere. Furthermore, the scalar-valued quantity ρ(*r*_{i},** e**)=ρ(

*r*_{i},−

**) characterizes a normalized odf that, in addition to the internal variables**

*e*

*r*_{i}, allows the material’s anisotropy to be modelled. For the general case, where , the evolution of this odf also contributes to the dissipation inequality.

## 4. A remodelling formulation

Remodelling, interpreted as a change in internal structure of the material, is reflected by evolution equations for the direction vectors *r*_{i}. Before we discuss a particular ansatz to incorporate these effects, additional assumptions on the anisotropic part of the Helmholtz free energy are addressed first. On the one hand, Ψ^{ani} is understood to be independent of the orientation (sign) of the direction vectors, in other words . Apart from this rather natural assumption, which similarly is well established for even-order structural tensors, we, on the other hand, directly relate the internal variables to the numerical framework, i.e. the integration of equation (3.4). As the algorithmic integration over the unit sphere reduces to a summation over a finite number of integration directions, say for , we model these unit-vectors to coincide with the internal variables so that *m*=*p* and *e*_{i}=*r*_{i} with . Based on this assumption, equation (3.4) is approximated as
4.1
with *w*_{i} denoting the integration factors, which apparently depend on the integration scheme chosen, and follows by analogy with equation (2.2).

In general, anisotropy evolution—in other words, remodelling—can be incorporated by means of evolving odf quantities. As this contribution proceeds, we combine this idea with making use of the well-established concept of evolving internal variables. In view of the particular microsphere model at hand, one may either consider the odf factors *w*_{i}ρ_{i} or the corresponding integration directions *r*_{i} as internal variables. Combinations of both approaches are possible as well—in this contribution, however, we assume *w*_{i}ρ_{i}=*w*_{i}=const.∀*r*_{i} and *r*_{i}≠const. To later on visualize local anisotropic material properties and as a macroscopic measure of anisotropy, a second-order generalized structural tensor is introduced
4.2
where *n*_{k}⋅*n*_{l}=δ_{kl} and *A*_{1}≥*A*_{2}≥*A*_{3}≥0 with . Even though initial anisotropy can generally be included, we will consider an initially isotropic setting so that ** A**|

_{t0}=[1/3]

**, whereby**

*I***is the second-order identity tensor. Note that deformation-induced anisotropy evolution requires in order to be able to recapture a stress-free reference state.**

*I*### (a) Constitutive model for the evolution equations

The quantity driving the remodelling process, as captured by deformation-induced anisotropy evolution, is here referred to as the reduced dissipation relation (3.1). In view of an individual integration direction—compare equation (4.1)—one obtains
4.3
The evolution of *r*_{i} is now motivated by the projection of this, say, driving force onto the plane perpendicular to *r*_{i} itself. Accordingly, the constraint for the rate of the unit vector *r*_{i} is automatically satisfied. Conceptually speaking, a shear-type deformation measure is considered by means of the projection operator ** I**−

*r*_{i}⊗

*r*_{i}. In this regard, we assume 4.4 wherein is a weighting factor. Note that is identically zero for

*r*_{i}being co-linear with one of the principal directions of . Based on this ansatz, the reduced dissipation inequality results in 4.5 which constitutes a semidefinite quadratic form that, by analogy with , identically vanishes in case

*r*_{i}corresponds to a principal direction of . As a typical example, we have for , while yields

*r*_{i}to align with the dominant principal direction of . The remaining task consists in particularizing

*f*

_{i}. Apparently, the standard form of the dissipation inequality for the closed solid mechanical system is satisfied for with

*c*≤0. For the problem at hand, however, a sound description of the adaptation of the biological tissue to be modelled is described by the theory of open systems, which is not outlined in this contribution—in addition, a chemical potential might, at least formally, be included as well. In this context, consider a sufficiently loaded body undergoing non-zero Dirichlet boundary conditions fixed in time. An adaptation process should then render the tissue to increase its loading capacity so that the system has to be fed with energy. Obviously, the reduced form of the standard dissipation inequality, as represented by , is violated in this case, which would similarly be exhibited for healing processes.

As the adaptation of biological tissues is usually bounded by certain, say biological, limits, the evolution equation for *r*_{i} should account for so-called saturation effects. On the one hand, the evolution of *r*_{i} saturates for *r*_{i} aligning with a principal direction of . On the other hand, we here restrict the maximum degree of anisotropy by making use of the spectral decomposition of ** A**. To be specific,

*r*_{i}is assumed to evolve as long as the difference between the largest and smallest eigenvalue of

**, i.e.**

*A**A*

_{1}−

*A*

_{3}, remains smaller than a pre-defined limit value, denoted here as 0<

*A*

_{Δ}≤1. In addition, a relaxation-type time parameter,

*t**>0, is incorporated and, moreover, we set in case the related fibre stretch does not exceed a certain threshold, ; note that for . Alternatively, one could also set up a dead-zone-type remodelling equation, for instance, by setting

*f*

_{i}=0 for taking values within a specific interval. In summary, the proportionality factor

*f*

_{i}introduced in equation (4.4) is assumed to be 4.6

### (b) Numerical integration of the evolution equation

In order to integrate the individual direction vectors *r*_{i} in time, we consider a finite time interval, Δ*t*=*t*_{n+1}−*t*_{n}≥0, and apply a Euler-type method, i.e.
4.7
together with to satisfy the normalization constraint. Alternatively, for example, rotation-parameter-based or exponential integration schemes could be adopted, which automatically guarantee .

Note that the remodelling formulation itself, respectively the evolution of deformation-induced anisotropy, depends on the integration scheme chosen for the unit sphere. Even though a perfect alignment of all integration directions with respect to one axis of anisotropy, say , renders ** A**=

**⊗**

*a***to be independent of the particular integration scheme, any intermediate state obviously depends on the type of the particular algorithmic formulation. Moreover, the orientation of the base system, with respect to which the fixed coordinates of the initial integration directions**

*a*

*r*_{i}|

_{t0}are referred to, influence the simulation results as well—in particular as direction vectors co-linear to a principal direction of , or in the respective plane spanned by two principal directions in case of two identical principal values of , do not reorient. When considering states of deformation with changing principal directions of deformation, however, this effect is almost negligible and, moreover, higher order integration schemes show common convergence properties.

Further issues of implementation—with a special focus on iterative finite element approaches—are discussed in appendix A.

## 5. Numerical examples

Before discussing some representative numerical examples, we first specify the particular form of the Helmholtz free energy. Based on this, homogeneous as well as inhomogeneous states of deformation are discussed. Additional details relevant for the finite element implementation to compute general boundary value problems are provided in appendix A.

### (a) Specification of the Helmholtz free energy

As commonly observed for soft biological tissues, the high water content of the material’s ground substance results in an almost incompressible response. In this context, we make use of a related quasi-incompressible neo-Hooke-type model for the volumetric and isotropic isochoric part of the Helmholtz free energy, namely
5.1
The anisotropic fibre-related part of the Helmholtz free energy is assumed to be characterized by means of a worm-like chain model. While the statistical and entropy-based background of this one-dimensional nonlinear model is extensively discussed in the references cited in § 1, we here rather adopt this approach. In addition, a repulsive term is incorporated in order to be able to generally recapture a stress-free contribution upon unloading, i.e. . In contrast to the neo-Hooke model applied for the isotropic part, the anisotropic fibre contributions are assumed to undergo solely affine deformations. Accordingly, the macroscopic stretch along *r*_{i} fully determines the anisotropic stress contribution in this direction. In summary, the worm-like chain model considered takes the representation
5.2
wherein the Boltzmann constant is denoted by *K*=1.38×10^{−23} (JK^{−1}), θ is the absolute temperature and *A* is established as the persistent contour length. While *r*_{0} characterizes the length of the chain for the undeformed state, i.e. at , the actual chain length follows from . From the mechanical point of view, it is worth mentioning that both the Helmholtz free energy and the related stress contribution both approach infinity for , which reflects one of the main ideas of the worm-like chain model. Furthermore, any response under compression has been excluded and, formally, the constant is incorporated so that *n* ψ=0 for . Note that one alternative could also have assumed the worm-like fibre contributions to depend on the total stretch, ** C**, rather than on the isochoric part, .

In view of the set of material parameters chosen, we adopt constants representative for blood vessels as summarized in Alastrué *et al.* (2009*a*). The particular parameters used are summarized in table 1, wherein the abbreviating notation *B*=*n**K*θ*r*_{0}/[4*A*] has been applied—the influence of varying material parameters is studied later on. Moreover, we will also use a Cartesian base system as characterized by the set of orthonormal vectors {*e*_{1},*e*_{2},*e*_{3}}, respectively {*e*_{x},*e*_{y},*e*_{z}}.

Concerning the numerical integration on the unit sphere, a simple algorithm is adopted based on 21 integration points for half of the sphere, i.e. *m*=42. A visualization of the underlying integration directions, together with a graphical illustration of the method of stereographic projection, is highlighted in figure 1. Detailed background information on this particular integration scheme, which nowadays can be considered as the standard approach for isotropic materials, is provided in the contributions by Bažant & Oh (1986) and Miehe *et al.* (2004). In view of higher order methods and modifications of standard integration algorithms, the reader is referred to, for instance, Heo & Xu (2001) and Alastrué *et al.* (2009*a*,*b*). As mentioned in §4, the numerical results strongly depend on the particular integration scheme used. In fact, it turns out that using *m*=42 cannot capture high degrees of anisotropy if *r*_{i}=const. However, the remodelling formulation proposed is less sentitive to the particular integration algorithm. As indicated above, an illustrative example consists of the fact that all suitable integration schemes render an identical distribution for the fully aligned case, i.e. ** A**=

**⊗**

*a***.**

*a*Apart from applying the method of stereographic projection also to other quantities of interest, such as stretches and stresses, we additionally make use of an odf-type representation. To be specific, the visualization of ρ^{A}=** e**⋅

**⋅**

*A***, for , allows visualization of the anisotropic material properties captured by the symmetric second-order tensor**

*e***.**

*A*### (b) Homogeneous deformations

In the following, the model proposed is investigated for different homogeneous deformations, namely isochoric equibiaxial tension, pure shear and simple shear. Special emphasis is thereby placed on the evolution of deformation-induced anisotropy, which we will illustrate by means of ** A** in terms of the odf-type function ρ

^{A}and via the difference between the maximal and minimal principal value,

*A*

_{1}−

*A*

_{3}. The particular loading history considered is based on linearly increasing the representative loading parameter, here either λ or γ, within a time period of 10 time steps and then fixing its value for a time period of 70 steps (figure 2). What we refer to as one step is related to one time unit—i.e. 1[s], which could be in the range of minutes, hours, days or weeks—so that the loading is increased within

*n*

_{0}=10 steps. However, we will additionally study rate dependencies during the loading process as the quasi-static model is of viscoelastic type. Accordingly, setting

*n*

_{0}=5 reflects twice as fast loading process as characterized by

*n*

_{0}=10, while

*n*

_{0}=20 reduces the speed of loading by a factor of two compared with

*n*

_{0}=10.

#### (i) Isochoric equibiaxial deformation

A state in isochoric equibiaxial deformation can be represented via the deformation gradient ** F**=λ

*e*_{1}⊗

*e*_{1}+λ

*e*_{2}⊗

*e*_{2}+λ

^{−2}

*e*_{3}⊗

*e*_{3}. Figure 3 shows the saturation-type behaviour of the deformation-induced anisotropy evolution by means of visualizing the difference between the maximal and minimal principal value of

**, i.e.**

*A**A*

_{1}−

*A*

_{3}. Moreover, the influence of individual material parameters is displayed by choosing values respectively different from those given in table 1.

In this regard, figure 3*a* displays the influence of the saturation value *A*_{Δ}. It is clearly seen that a value of *A*_{Δ}>1 has almost no influence when compared with *A*_{Δ}=1—the value to *A*_{1}−*A*_{3} may only slightly increase faster. Nevertheless, *A*_{Δ}=1 constitutes an upper limit for *A*_{1}−*A*_{3} so that also large values of *A*_{Δ} render the same saturation plateau. It is also worth mentioning that the extremal value, i.e. *A*_{1}−*A*_{3}=1, is not approached here for three reasons: (i) the equibiaxial deformation can (for ) at most result in *A*_{1}−*A*_{3}=1/2, (ii) directions undergoing a stretch smaller than do not reorient, and (iii) the integration directions co-linear with the principal stretch direction *e*_{3} remain constant. In view of the value *A*_{Δ}=0.2, one particularly observes that the anisotropy evolution takes place in a viscous manner as the loading parameter λ is fixed after 10 steps and the graph of *A*_{1}−*A*_{3} continues to increase.

The influence of the parameter is highlighted in figure 3*b*. As expected, setting ends up with the largest saturation value for *A*_{1}−*A*_{3} as also directions undergoing all levels of compression reorient. In case reorientation under compression is excluded, however, the different values of here result in the same saturation value for *A*_{1}−*A*_{3} as the integration scheme applied to the unit sphere is rather coarse.

Since the remodelling formulation proposed is time dependent, namely of viscoelasticity type, the time interval within which the respective loading is applied strongly influences the adaptation process. By analogy with strain rate dependencies, figure 3*c* displays the influence of the time period chosen for the loading interval. Apparently, a higher, say, strain rate yields the remodelling process to develop faster. Similarly, one could also modify the relaxation-time-type parameter *t** (figure 3*d*). Obviously, larger values of *t** reduce the remodelling velocity.

The deformation-induced anisotropy evolution is additionally pictured in figure 4, in which an odf-type visualization for the symmetric second-order tensor ** A** has been applied together with the method of stereographic projection. On the one hand, we observe for the three different states of deformation (figure 4

*a*–

*c*) that the degree of anisotropy continuously evolves with time as the odf increasingly deviates from its initially spherical distribution. On the other hand, the equibiaxial state is clearly recovered by the, say, circular cross sections of the odf. To be specific, the equibiaxial loading results in a uniaxial distribution of ρ

^{A}so that

**possesses two identical principal values. To give an example, for step 80, we obtain**

*A**A*

_{1}=

*A*

_{2}=0.4034 and

*A*

_{3}=0.1933.

#### (ii) Pure shear

Similar to a state in equibiaxial tension, the principal stretch directions of a homogeneous deformation under pure shear, i.e. ** F**=λ

*e*_{1}⊗

*e*_{1}+λ

^{−1}

*e*_{2}⊗

*e*_{2}+

*e*_{3}⊗

*e*_{3}, do not depend on the deformation measure λ.

The saturation-type behaviour of the scalar value *A*_{1}−*A*_{3}, which represents the degree of deformation-induced anisotropy, is shown in figure 5. By analogy with figure 3, the influence of the different parameters is studied and, as expected, a similar response as exhibited for the equibiaxial deformation is observed. An odf-type representation is displayed in figure 6, where it is clearly shown that the degree of anisotropy increases with time. Compared with the previous example in §5*b*(i)—in which a uniaxial distribution was obtained with the two identical principal values of ** A** being larger than the remaining one—the distribution of ρ

^{A}here turns out to be biaxial. The second-order tensor

**, to which ρ**

*A*^{A}is referred to, possesses three different principal values; to be specific, step 80 renders

*A*

_{1}=0.6392,

*A*

_{2}=0.2332 and

*A*

_{3}=0.1276. Note, however, that a perfect alignment for and

*r*_{i}≠±

*e*_{3}∀

*i*results in a uniaxial distribution with

**=**

*A*

*n*_{1}⊗

*n*_{1}.

#### (iii) Simple shear

The elaboration of a homogeneous deformation in simple shear, i.e. ** F**=

**+γ**

*I*

*e*_{1}⊗

*e*_{2}, is of special interest for investigations on anisotropic constitutive equations since the principal stretch directions depend on the deformation measure γ.

By analogy with the two homogeneous deformations previously discussed, figure 7 highlights the saturation-type behaviour of *A*_{1}−*A*_{3}, whereby similar effects as before are recovered. Furthermore, the odf-type representations in figure 8 show that the integration directions align according to the dominant principal stretch direction. The distribution itself reflects that ** A** possesses, in general, three different principal values so that ρ

^{A}is of biaxial form. In line with the homogeneous deformation in pure shear, the principal values of

**for step 80 are**

*A**A*

_{1}=0.5748,

*A*

_{2}=0.2512 and

*A*

_{3}=0.1740.

### (c) Inhomogeneous deformation

Motivated by the experimental investigation reported in Eastwood *et al.* (1998), we next qualitatively study the inhomogeneous deformation of a plate loaded under tension. The particular specimen considered possesses the dimensions 80×40×4 and its finite element discretization has been performed with 20×10×2 enhanced eight node bricks, Q1E9, as advocated by Simo & Armero (1992). Concerning the set of material parameters, the values summarized in table 1 are adopted except for and *L*=2. Moreover, quasi-static loading conditions are applied by means of Dirichlet boundary conditions so that the plate is longitudinally stretched to a length of 120—the transverse displacements at these boundaries being clamped. By analogy with the homogeneous deformations discussed in § 5*b*, the respective load increments are linearly increased within 150 steps, whereas the boundary conditions thereafter remain fixed for a period of another 150 steps; related graphical illustrations are highlighted in figures 9 and 10. The sum of reaction forces present at the respective Dirichlet boundaries are given in figure 9. As expected, we observe that the reaction forces continuously increase. While the slope of this curve increases during the loading interval (steps 1–150), the amount of reaction forces tends to saturate for the interval within which the prescribed degrees of freedom are fixed (steps 151–300).

Figure 11 displays the distribution of the Cauchy stresses with respect to the macroscopic longitudinal loading direction for different loading levels. Apparently, the stress values increase with deformation even for the loading interval within which the prescribed displacements are fixed. Practically speaking, this effect cleary underlines the viscous anisotropy evolution. Please note that the individual contour plots are referred to by different legend scales. As experimentally observed, the fibre distribution aligns according to the local loading direction and the degree of anisotropy increases towards the middle necking zone. This effect is recaptured by the proposed remodelling formulation and is shown in figure 12, where contour plots of the anisotropy measure *A*_{1}−*A*_{3} at different steps of the simulation underline the inhomogeneity of the tension problem considered as well as the almost isotropic material behaviour in the region of the bottom and top surfaces. Finally, figure 13 illustrates the angle between the macroscopic longitudial loading direction and the principal direction *n*_{1} related to the dominant principal value *A*_{1} of the generalized structural tensor ** A**. As identified from experimental measurements, the proposed modelling framework reflects an almost perfect alignment of the fibres with respect to the dominant stretch direction.

Nevertheless, at this stage, the simulation provides mainly qualitative results as the calculation of precise quantitative results requires further information on the experimental data such as integral load–displacement curves, local measurements of fibre distributions and so forth.

## 6. Summary

Continuum modelling approaches suitable for the simulation of biological tissues should account for the material’s relevant micromechanical properties and, depending on the particular application of interest, the dominating biological effects of the body considered. As various soft biological tissues are composed of a ground substrate and fibre assemblies such as collagen, we here additively combined an isotropic stress term with stresses related to the contribution of fibres, which in general render the mechanical constitutive law to be anisotropic. The formulation of these fibre contributions was based on a worm-like chain model as generally derived from a sound statistical and micromechanically motivated background. By adopting a computational microsphere-based strategy, this one-dimensional formulation could be applied within a three-dimensional context. Apart from the purely passive response, the biological adaptation of the biological tissues was incorporated as well. Special emphasis has thereby been placed on the formulation of a remodelling framework while growth phenomena are not included in the model at this stage.

In this work, remodelling is understood to be a process that renders the internal substructure of the material to be able to adapt to the local loading conditions. For the particular model proposed, we, accordingly, assume the anisotropic part to be the stresses to account for these phenomena. In other words, the referential orientations of the individual fibre contributions depend on the local loading conditions and change in a viscoelastic manner with time. Such an alignment of fibres is often also denoted as reorientation or, from the biological point of view, as turnover. The model developed directly combines this remodelling with the computational microsphere approach. To be specific, the respective directions introduced to perform the numerical integration over the unit sphere are reoriented. As a result, the formulation accounts for deformation-induced anisotropy evolution. Without loss of generality, we adopted a distribution of integration directions that yields the initial body to be isotropic. In general, however, this distribution can be weighted by an initial odf so that initial anisotropy can be captured as well. Note that a similar algorithmic approach could also be applied to hard biological tissues.

The particular evolution equation, which described the reorientation of a single integration direction, was assumed to be driven by a representative stretch measure and rendered the related mechanical dissipation term to take a semidefinite quadratic form. Alternatively one could also choose a, for instance, stress-driven ansatz. Moreover, the evolution equation proposed results in an alignment of the integration directions with the dominant principal stretch direction so that the biological tissue maximizes its load-carrying capacity. Saturation effects are, on the one hand, naturally included as the remodelling stops as soon as a direction is aligned with a principal stretch direction. On the other hand, an additional saturation value has been introduced to be able to further limit the maximal degree of anisotropy of the tissue. The numerical examples investigated showed the basic algorithmic applicability of the modelling framework and captured fundamental phenomena such as reorientation and remodelling effects observed for soft biological tissues. The formulation developed can be applied to the simulation of general boundary value problems and, owing to the fact that the initial body has been assumed to be isotropic, can also be used to estimate or rather detect the material’s anisotropic properties under physiological loading conditions.

Concerning future research, the design of evolution equations that on the continuum level capture fibre alignment with more than one single direction is of special interest as several biological tissues show macroscopically orthotropic behaviour. Furthermore, the investigation of different driving forces and a possible combination with macroscopic target stresses constitute an attractive line of research. Apart from the isotropic stress contributions, the anisotropic terms were based on an affine microsphere model—in future, research on the influence of non-affine anisotropic stress terms is of interest as well. An issue of key importance is the combination of the remodelling formulation proposed with growth phenomena. In view of the algorithmic implementation, a particular integration scheme has been applied to the unit sphere, which turned out to be part of the model itself and directly influenced the simulation results. Future investigations should elaborate different integration schemes and focus on their efficient numerical treatment. Alternatively, direct orientation-distribution-based approaches can be adopted for the modelling of anisotropic biological tissues—for example, in combination with modified integration schemes for the unit sphere or in terms of an additional balanced equation.

## Acknowledgements

A.M. gratefully acknowledges the partial financial support by the Swedish Research Council (Vetenskapsradet) under the grant 2007-5224.

## Appendix A. Some aspects of implementation

The equilibrium of a single-phase solid continuum can be represented by the balance of linear momentum, which for the quasi-static case and in the absence of volume forces takes the local representation
A 1
Based on equations (3.3) and (4.1), the hyper-elastic form of the Piola–Kirchhoff stresses, ** S**=2 ∂

_{C}Ψ, results in A 2 Apparently, the detailed representation of the fourth-order isochoric projection operator thereby reads A 3 with I

^{sym}being the fourth-order symmetric identity tensor. In view of the micro-sphere approach adopted, compare equation (4.1), the anisotropic stress contribution is further specified as A 4 In order to abbreviate notation, the index •

^{n+1}has been neglected—a simplification we also apply in the following.

As a Newton-algorithm-based finite-element approach requires the computation of the corresponding tangent operator, we here mention its ‘material’ part, E=2d_{C}** S**, which takes the form
A 5
By analogy with equation (A 4), we next place emphasis on the anisotropic term of the tangent operator
A 6
While the computation of is straightforward, as well as and , we discuss the last two contributions of equation (A 6) in more detail. On the one hand, one observes
A 7
with for

**being a second-order tensor, whereas**

*T***and**

*v***are vectors. On the other hand, one has to additionally account for the update and normalization when calculating the remaining derivative. To be specific, equation (4.7) yields A 8 With these results in hand, the contribution of an individual integration direction to the anisotropic part of the tangent operator results in A 9 Owing to the non-associated character of the evolution equation (4.4), the tangent operator possesses only minor symmetry.**

*w*## Footnotes

One contribution of 12 to a Theme Issue ‘Mechanics in biology: cells and tissues’.

- © 2009 The Royal Society