## Abstract

The importance of tissue remodelling is widely accepted, but the mechanism by which the remodelling process occurs remains poorly understood. At the tissue scale, the concept of tensional homeostasis, in which there exists a target stress for a cell and remodelling functions to move the cell stress towards that target, is an important foundation for much theoretical work. We present here a theoretical model of a cell in parallel with a network to study what factors of the remodelling process help the cell move towards mechanical stability. The cell-network system was deformed and kept at constant stress. Remodelling was modelled by simulating strain-dependent degradation of collagen fibres and four different cases of collagen addition. The model did not lead to complete tensional homeostasis in the range of conditions studied, but it showed how different expressions for deposition and removal of collagen in a fibre network can interact to modulate the cell's ability to shield itself from an imposed stress by remodelling the surroundings. This study also showed how delicate the balance between deposition and removal rates is and how sensitive the remodelling process is to small changes in the remodelling rules.

## 1. Introduction

Human tissues routinely change in response to environmental cues. The fundamental processes of growth and remodelling are necessary for the development of healthy tissues and for healing of tissues in the aftermath of injury or illness. Growth and remodelling can occur by a wide range of mechanisms, including cell proliferation [1] or death [2,3], pattern formation [4,5] and synthesis and/or degradation of extracellular matrix (ECM) [6–8]. This last process, cell-driven modification of the ECM, is of great interest because of the critical mechanical role played by ECM proteins, especially collagen, in many tissues [9]. The collagen network is often responsible for maintaining mechanical integrity and tensile strength of the tissue [7,10–15]. Although many other components (e.g. elastin [16], proteoglycans [17]) are also essential for tissue function, the ubiquity of collagens in load-bearing tissues makes them an attractive target for theoretical and experimental study. Collagen remodelling reflects the complex biomechanical relationship between the ECM and the cells that it supports.

In order to understand the remodelling of collagenous tissues, the underlying biochemical and biomechanical interactions between cells and the ECM must be explored. Significant strides have been made through the use of theoretical models that treat the cells and ECM as a continuum (e.g. [18–22]). In particular, many models invoke the concept of tensional homeostasis, in which cells grow and/or modify the ECM to achieve a target stress. In its simplest form, the tensional homeostasis principle postulates that a target stress *σ*_{o} exists, and that some remodelling process occurs at a rate proportional to the deviation (*σ* − *σ*_{o}) between the cell stress and that of the target. Under most circumstances, the target stress is stable—actions taken in response to a deviation tend to reduce that deviation and return the system to the homeostatic stress state. Examples of this principle in models of tissue remodelling may be found throughout the biomechanics literature [23–26], and it is obviously related to the effect of matrix mechanics on cell behaviour [27–29].

It remains unclear, however, how tensional homeostasis is achieved or what factors are necessary for tensional homeostasis. Tissue-level stress must be translated into tensions on the individual cells, and any tissue-level response must be driven by individual cellular responses. The details of these translations and the mechanisms by which the cell changes the collagen fibres around it are far from obvious, and the complex interactions between the different individual components of the tissue make exploration of these questions difficult. Experimental studies have shown, for example, that collagen synthesis is altered by the stretch [30] and the stiffness [31] of the substrate on/within which fibroblasts are cultured, and the enzymatic degradation of collagen is dependent on the collagen fibre tension [32]. The question at hand, then, is whether remodelling of the ECM through synthesis and degradation of collagen fibres can lead a system to tensional homeostasis. We found that proportional remodelling alone did not lead to a fully stable system, but that it was an efficient method to move towards stability with minimal feedback required from the cell.

## 2. Material and methods

To study the stresses on a remodelling system of collagen networks and cells, we used an idealized model of a network and cell in parallel (figure 1). In this configuration, the cell and network experience the same strain, and their individual stresses sum to the total stress of the system. The system was stretched by varied amounts and held at constant stress during remodelling. While the schematic in figure 1 shows a single cell, the ‘Cell’ component of the system is representative of the collection of cells in a tissue.

Figure 2 presents an overview of the multistep remodelling strategy. First, the fibre network was generated. Next, the system, consisting of a neo-Hookean cell and a fibre network in parallel, was stretched to a specific stretch level, and the system stress was recorded. The system was then allowed to remodel by adding and/or removing material, and the strain was adjusted to maintain the same applied stress. The remodelling and adjusting steps were repeated until the entire system had equilibrated, and the final cell and network stresses and strains were recorded. Details on each step in the process are provided in the subsequent sections.

### 2.1. Collagen network and cell models

Collagen networks were modelled as three-dimensional networks of nonlinear hyper elastic springs connected at freely rotating nodes as we have done previously [33]. The networks were generated by randomly seeding points inside a unit cubic representative volume element (RVE), and Delaunay triangulation was used to connect polygonal elements from the seed points. The edges of the polygons became the collagen springs and the vertices became the nodes. Five random Delaunay networks were created with 714 ± 48 (mean ± s.d.) collagen fibres. The fibres were assigned an initial radius of 50 nm and stiffness of 79 MPa [34]. The fibre volume fraction of the RVEs was 1.2 ± 0.07%. The real size of the RVE was determined by a conversion factor between the computational and real domains that is defined in terms of the fibre radius. Setting the collagen fibre radius to 0.005 computational units (CUs) gave us a desired collagen volume fraction of approximately 1%. Using the conversion 0.005 CU = 50 nm makes the real size of the RVE 10 µm.

The cell was modelled as an incompressible neo-Hookean solid which was enforced by setting the Jacobian J = 1, preserving the volume. The Cauchy stress on the cell, was then defined as
2.1
where *C* is a constant defining the stiffness of the cell. Depending on the cell type, cell stiffness can vary widely; therefore, the value of *C* was defined with respect to the stiffness of the unremodelled collagen network to explore three different cases:

(i) a cell twice as stiff as the network—could also represent a system of low-density matrix and many cells;

(ii) a cell with the same stiffness as the network—cell and matrix have the same contribution in the system; and

(iii) a cell half as stiff as the network—could also represent a system of high-density matrix and few cells.

Network stiffness was defined as the slope of a secant line on the stress–strain curve from the origin to half-maximum stretch (before remodelling).

### 2.2. System deformation

The cell-network system was stretched by different amounts (5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55 and 60% strain) before remodelling. The system was stretched uniaxially and incompressibility was enforced. Stretching was performed as previously described [35]. Briefly, the boundaries of the network were moved in steps of 1% in the direction of stretch, which we will call the 1-direction. In the other two directions, the stretch ratio was set to preserve the volume of the RVE (*λ*_{2} = *λ*_{3} = 1/√*λ*_{1}). The internal network nodes were iteratively moved until all were at mechanical equilibrium. The force, *F*, on each node was defined using an exponential relation such that the forces on fibres in tension are much larger than forces on fibres in compression with a computationally advantageous continuous force function.
2.2
where *E*_{fibre} is the collagen fibre modulus at infinitesimal strain, *A*_{fibre} is the cross-sectional area of the fibre, *ɛ*_{G} is the Green strain of the fibre and *B* is a nonlinearity parameter which was set to 1.2 based on previously fitted data from collagen gels [36].

The volume-averaged stress on the equilibrated network, was calculated as
2.3
where is the volume-averaged Cauchy stress tensor, *V* is the current volume of the RVE, *bc* indicates discrete summation over the boundary nodes, *x _{i}* is the position of a node on the boundary,

*f*is the force on that node,

_{j}*p*is the hydrostatic pressure term for incompressible materials and

*δ*is the identity. The value of

_{ij}*p*was obtained by solving for the stresses on the free surfaces. For full derivation, refer to [37].

### 2.3. Network remodelling

Collagen network remodelling was represented in the model by increasing or decreasing the radius of the cylindrical collagen fibres based on addition and removal rates, respectively. The radius-based approach does not allow creation of new fibres or the complete removal of a fibre.

#### 2.3.1. Collagen removal

The rate of fibre removal, was related to the stretch of the fibres such that more highly stretched fibres were more protected from degradation, following the observations of [32]:
2.4
where *λ*_{fibre} is the stretch ratio of the fibre, *r* is the fibre radius and *r*_{small} is a value representing a small radius, chosen here to be 10% of the original value. The factor ensures that the fibre radius will not become negative during the course of remodelling. When the ratio approaches unity and (d*r*/d*t*) becomes *r*-independent, leading to linear decrease in the case of constant stretch. When however, the ratio approaches *r*/*r*_{small}, and the radius decays exponentially to zero. The scaling constant *k*_{1} has units of nanometre per second, and *k*_{2} is a dimensionless constant determining the amount of exponential dependence on stretch. The negative exponential form of the governing equation gives more stretched fibres (larger *λ*_{fibre}) lower degradation rates.

The values of the constants *k*_{1} and *k*_{2} were obtained previously by fitting the experimental data of [32]. Hadi *et al*. [34] found *k*_{1} to be 1.05 × 10^{−3} nm s^{−1} and *k*_{2} to be 0.83. The obtained values are for a specific *in vitro* case of collagen degradation by bacterial collagenase. To look at a more general case of collagen remodelling, we write the equation in dimensionless form. If we let the characteristic time scale to be *τ* = *R*/*k*_{1}, where *R* is the initial value of the fibre radius (50 nm), then the scaling factor for the removal rate can be re-written as follows:
If we further define the dimensionless constants and then the dimensionless rate of removal becomes
2.5

This degradation expression was used in all simulations.

#### 2.3.2. Collagen addition

A vast array of options exists for modelling fibre deposition. In this paper, we begin with a general form based on the principle that two factors can influence deposition: (i) the cell may tend to deposit more material when it is more stretched (note that because of the elastic constitutive model for the cell, this assumption is equivalent to a cell stress-based model) and/or (ii) the cell may tend to align and thus deposit new matrix in the direction in which the current fibres are already aligned (via contact guidance [38–40]). Assuming separability of the two factors, we write the general form
2.7
where *k*_{3} is a constant.

For the cell stretch effect, we postulate an exponential factor such that
2.8
where *λ*_{cell} is the cell stretch ratio and is equal to the matrix stretch under the assumptions of figure 1. The constant *k*_{4} captures the sensitivity to stretch. A large value of *k*_{4} indicates a very stretch-sensitive cell, and *k*_{4} = 0 indicates stretch-independent deposition. For the fibre orientation effect, we scaled the deposition rate according to a von Mises distribution
2.9
where *θ* is the orientation of a fibre to be deposited, *θ*_{o} is the average orientation of the existing matrix, *I*_{0}(*κ*) is the modified Bessel function of order zero, and *κ* is a shape parameter that indicates stronger alignment for larger *κ*. If *κ* = 0, *g*(*θ*) = 1 for all *θ*, and all directions have equal strength; as *κ* goes to infinity, *g* approaches a delta function at *θ* = *θ*_{o}. Because of the uniaxial stretch in this study, the direction of *θ* was always the 1-direction. The normalization factor of *π* normally present in the von Mises distribution has been adsorbed into the constant *k*_{3}. This model does not account for the degree of orientation in the matrix, only the direction; a strongly aligned network is no more likely to induce aligned deposition than a weakly aligned one.

Combining these expressions, the total deposition rate is given by
2.10
or, in the dimensionless form following equation (2.5):
2.11
where *α* = *k*_{3}/*k*_{1} is a dimensionless quantity representing the rate of deposition relative to removal.

In the current study, we considered four cases (table 1):

(1) Stretch-independent, isotropic deposition (

*k*_{4}=*κ*= 0).(2) Stretch-dependent, isotropic deposition (

*k*_{4}≠ 0,*κ*= 0).(3) Stretch-independent, anisotropic deposition (

*k*_{4}= 0,*κ*≠ 0).(4) Stretch-dependent, anisotropic deposition (

*k*_{4}≠ 0,*κ*≠ 0).

The dimensionless parameters were varied to cover a broad range of possibilities. The values used in each case are presented in table 1.

The net rate change of radius based on the above expressions was 2.12 and the new fibre radius was calculated by a forward Euler step, 2.13

### 2.4. Remodelling at constant stress

When the radius for each fibre was updated, the stress in the system changed. To return the system to the applied stress, the strain was adjusted (either by contracting the system or by stretching it further). The adjustment was done by iterating the strain until the new stress was within a specified tolerance of the applied stress (±1 × 10^{−5} kPa). Once the applied stress had been achieved, another remodelling time step was taken, and the rates of removal and addition were calculated for each fibre again. A total of 3000 remodelling iterations were performed (each iteration representing 20 s remodelling time), enough for the cell and network stresses to stabilize. The equilibrium cell and total system stresses were plotted versus each other as shown in figure 3.

## 3. Results

The deformation steps before remodelling starts were the same for all cases. Cellular and network stiffness varied in the range of 150–300 kPa before remodelling. For the case of a cell with half the stiffness of the network, the cell stiffness range was between 75 and 150 kPa. For the stiff cell cases, cell stiffness was over 300 kPa.

### 3.1. Case 1 (*k*_{4} = *κ* = 0)

The rate of collagen removal in this case, as in all others, was based on the stretch of the fibres. The rate of addition was set to the dimensionless constant *α*. The value of *α* had an upper limit of 1, for which the rate of addition was equal to the rate of removal at zero stretch, which caused the network to deposit collagen continuously and prevented it from equilibrating (figure 4). On the other hand, a low *α* value of 0.8 led to more removal than addition. Net loss of collagen is reasonable, as would be expected in a tensional homeostasis model if the cell stress were below the target stress, but for this study, we considered systems that would lead to net collagen deposition. Therefore, the values of *α* were selected to be 0.9, 0.92 and 0.95.

At *α* = 0.9 and a medium stretch of 20%, the system maintained a constant stress of 171 kPa during remodelling as the cell and network stresses moved in opposing directions until both reached equilibrium (figure 5). The stress on the network increased during remodelling, prompted by an increase in the collagen volume fraction. The stress on the cell, in complement, decreased as the system stretch ratio decreased to maintain a constant total system stress. The volume fraction increase was a result of the net deposition of new material. The mean fibre radius in the network increased—all fibres started with a radius of 50 nm and during remodelling some fibre radii became much bigger (figure 6). The results presented with these parameters are representative of the qualitative behaviour seen for other values of *α* and for other stretches. The value of *α* had a significant impact on the equilibrium stresses (figure 7). An increase in *α* caused an increase in the collagen volume fraction, as expected, and it also corresponded to a decrease in the fraction of the total stress that was borne by the cell. For all three values of *α*, however, the cells never reached a constant stress value independent of the total stress. For example, in the low-slope (high-stress) region of the plot in figure 7*b*, a 16.6% increase in total stress led to a 5.4% increase in cell stress when *α* = 0.9. When *α* = 0.95, increasing the total stress by the same amount increased the cell equilibrium stress by 5.5%.

To check whether tensional homeostasis can be achieved by a cell of different stiffness, we varied the stiffness of the cell in relation to the collagen network, making it twice and half as stiff as the network (figure 8). The stiffer cell saw an increase in cell and total stress, as well as the fraction of the stress carried by the cell. With the softer cell, the cell was responsible for a smaller fraction of the total stress. Still, in all cases of cell stiffness, the cell stress became less sensitive to total stress but did not reach a clear plateau. At the higher stress, a 17.5% increase in total stress led to a 5.3% increase in cell stress for a soft cell. For a stiff cell, a 15.8% increase in total stress led to a 4.7% increase in cell stress.

### 3.2. Case 2 (*k*_{4} ≠ 0, *κ* = 0)

In Case 2, the cell was given more control of the remodelling process by allowing it to deposit more collagen the more stretched it became. Similar to Case 1, the cell stress decreased with remodelling while the network stress increased (figure 9). The plots in figure 9 are again for *α* = 0.9, but now if we look at two different system strains, 5% and 35%, we see that the stresses equilibrate much faster at higher strain. The larger the initial system deformation, the faster the cell responded because of the cell stretch-dependent deposition term. The cell bore only about 5% of total system stress. Compared to Case 1, the final cell stress was decreased in half.

### 3.3. Case 3 (*k*_{4} = 0, *κ* ≠ 0)

In Case 3, the cell was allowed to remodel the network by preferentially depositing collagen fibres in the direction of stretch but without changing the total deposition rate with stretch. The level of anisotropic deposition was controlled by the parameter *κ*. The values selected for *κ* were 0.02, 0.04, 0.06 and 0.08. The values were chosen again to ensure that the system would reach equilibrium, as larger values for *κ* caused a dramatic increase in the collagen volume fraction.

The representative example in figure 10*a* is with *α* = 0.9, *κ* = 0.2 and 20% initial strain. Similar to Cases 1 and 2, the cell stress decreased and network stress increased while the total system stress remained constant. Increasing *κ* increased the collagen volume fraction and decreased the system stretch ratio, similar to the effect of *α*. The fraction of the cell stress out of total system stress also decreased with increasing *κ* but the cell stress did not plateau for any of the *κ* values. For *κ* = 0.02, at high stress, changing the total stress by 16.6%, increased the equilibrium cell stress by 5.3%. For *κ* = 0.08, the cell stress increased by 4.8% for the same total stress increase as above. How the new fibres were deposited can be seen in figure 11. The fibres oriented in the 1-direction have a larger radius, though the effect of the increasing *κ* is not immediately obvious.

### 3.4. Case 4 (*k*_{4} ≠ 0, *κ* ≠ 0)

In Case 4, the cell was given the most control over the remodelling process, allowing it to deposit more collagen the more stretched it became and at the same time to deposit fibres in a preferred direction. The representative plot in figure 12 is for *α* = 0.9, *κ* = 0.02 and 20% initial stretch where the system equilibrates quickly. Increasing *κ* led to a small increase in the collagen volume fraction and a small decrease in the system stretch ratio. With the larger *κ* = 0.08, the cell had a smaller fraction of the total stress, compared to smaller *κ* = 0.02. The value of deposition at preferred orientation had a smaller effect on the remodelling process than the initial stretch (figure 13). The collagen fibre radius visibly increased when the initial stretch was changed from 5 to 35%, but the differences between Cases 2 and 4 were harder to see. Overall, Case 4 had the smallest fraction of the cell stress out of the total system stress (figure 14). At high stresses, a 16.6% increase in total stress led to a 4.6% increase in cell stress, when *κ* = 0.02.

## 4. Discussion

Tensional homeostasis has been a subject of interest for many years now, though experimental limitations have left many open questions [23]. The theoretical model presented here serves as a tool to explore how a cell can control its environment by remodelling its ECM. To model turnover of collagen fibre networks, we modelled the rate of collagen removal as dependent on the stretch of the fibres and we modelled the rate of addition in four different ways: (i) constant rate of addition, (ii) a more stretched cell deposits more collagen material, (iii) collagen material is deposited in the direction of stretch and (iv) the cell is sensitive to both stretch and fibre orientation (combination of Cases 2 and 3). The model allowed us to investigate how a cell might control its mechanical environment in an effort to achieve tensional homeostasis.

For all remodelling cases, the equilibrium cell stress approached a plateau but it never fully reached a zero slope. In general, increasing the rate of collagen addition led to a decrease in the fraction of the total system stress carried by the cell and an increase in the contraction of the system, moving it closer to its original configuration. The system was very sensitive to even small changes in *α* and could be induced into entering a non-converging state for larger values of *α*. Changing the stiffness of the cell in relation to the network changed the amount of cellular contribution to the system stress—a stiffer cell held a higher fraction of the total stress than a softer cell. The system was also very sensitive to changes in the initial stretch especially in cases where the cell could sense stretch and adjust its collagen deposition accordingly. The values used for *κ* in Cases 3 and 4 had to be kept very low, because, as with large *α*'s, the system did not reach a stable equilibrium for larger *κ*'s. An increase in *κ* led to an increase in the amount of added collagen and the contraction of the system. Even for very small *κ*'s, the system showed some improvement in reducing the cell stress fraction from the total stress compared to Case 1. The combination of allowing the cell to both sense stretch and align the collagen matrix resulted in the biggest reduction to cell stress compared to Case 1.

The trends observed when *α* and *κ* were large are best explained by recognizing that there is a minimum deposition rate and a maximum removal rate. For the isotropic cases (1 and 2), the dimensionless deposition rate must be at least *α*, and for the anisotropic cases, the dimensionless deposition rate for fibres oriented in the 1-direction goes as *α*(1 + *κ*) for small *κ* values. Because the sample is in tension, the minimum fibre stretch for most of the fibres is 1, and thus the minimum dimensionless removal rate is 1. Thus, if *α* > 1 or if *α*(1 + *κ*) > 1, then the deposition rate will *always* exceed the removal rate, and the fibre will thicken without bound.

Based on figure 14, it appears that very different modes of operation for the addition process can all be adequate for bringing the cell close to tensional homeostasis. The differences between Case 1and Case 4 were not great in terms of the final slope of the curves, but they were very different in the value of the cell stress—in Case 4, the cell stress was three times smaller. When the cell had more control over its remodelling environment, it was not only able to further reduce the stresses acting on it, but it was also able to reach remodelling equilibrium quicker. These factors could be important in protecting the cell from entering a pathological state, such as stiffening of the cell (as seen in tumour cells) [29].

There are certain assumptions and simplifications in our model that could merit further discussion. The cell was allowed to change only the collagenous matrix and not its internal structure; that is, the cell component was modelled as a homogeneous incompressible neo-Hookean solid with constant properties, whereas it is known that cell stiffness can in fact change in response to many factors, including the substrate stiffness [41,42]. Furthermore, in the current model, the cell was treated as entirely passive, with no active (contractile) contribution to the mechanical environment. New collagen fibres were deposited pre-stressed at a level equal to that of the already existing fibres. Finally, the remodelling feedback controls were based on a single proportional feedback rule. All of these assumptions are examined further below.

Our model did not take into consideration the internal structure of the cell, nor accounted for the many changes that can happen to it during remodelling. The exact molecular pathways of how mechanical stresses and strains are sensed by the cell and translated into a response are still not well understood. In this model, the cell was considered only as an elastic solid at three different stiffness levels. It is known that cells can change their stiffness with changes in stretch by reorganizing cytoskeleton components [43]. It is possible that giving the cell another level of control would change the equilibrium stresses. Mizutani *et al*. [43] found experimentally that cells returned to their pre-deformation stress levels in just 2 h, indicating that the cells can maintain an equilibrium stress level without the remodelling of the ECM. However, remodelling of the cytoskeleton and the ECM seem to occur on different time scales—with internal cellular remodelling occurring very rapidly compared to slower turnover of collagen fibres.

Of course, remodelling is a very complex process involving a multitude of chemical and mechanical pathways. Cells can release growth factors, hormones and enzymes, produce matrix and change cell alignment. The degradation of the collagen matrix can also be regulated by a number of factors, not just strain-dependent protection. Our model focused on the mechanical aspects of remodelling and how it applies to the local microenvironment of a cell. The different types of pathways available to a cell give it the ability to respond to different types of stimuli. Though mechanical and chemical stimuli cannot be fully separated, it is still informative to look at the mechanical response of a cell to mechanical cues. The mechanical responses that we covered have been appreciated rather recently and also contain some unanswered questions, for example how new collagens are integrated into the matrix. There is yet no experimental evidence on the stress state of newly deposited collagen, but for the cells to maintain a constant stress during remodelling, collagens would need to be deposited pre-stressed. In our work, new collagens were added by increasing the radius of already existing collagen fibres, which inherently meant that new collagens were pre-stressed.

The lack of complete tensional homeostasis in any of the models is reminiscent of the offset that exists in a proportional control system. Although the system in the current study is considerably more complex than a simple linear controller, the similarities are informative. The deposition rate is constant, so equilibrium for a given load must occur when the fibre stretch is such that the deposition rate matches the removal rate. This stretch, however, is an equilibrium stretch, not the set point stretch (i.e. the stretch at which no control action is taken, which in the current model would be zero stretch). The equilibrium is offset from the set point just as in a classical proportional control problem, and, as in the classical case, the offset is a function of the operating conditions, which in our system means that the equilibrium cell stretch must be a function of the imposed stress. A set point could be achieved through an adaptive control algorithm or through the addition of integral control, which may well exist for this purpose inside the cell and which would allow the deposition rate to depend not only on the stress but on the accumulation of some stress marker over time. The possible role of integral control within the cell is a topic of considerable study [44–47], and it would be interesting to explore it in the context of remodelling.

In conclusion, the results presented in this work show that, at least in some cases, a leveling-off similar to tensional homeostasis can arise naturally from the balance of deposition and removal without requiring that the cell ‘knows’ what the target stress is. Most importantly as relates to this issue, in our Case 1, a tensional homeostasis-like response (defined as small changes in cell stress over a wide range of tissue stresses because of remodelling) was achieved without any presumed sensory mechanisms in the cell. Even in the more complex models of our other cases, the remodelling was not driven by a desire to achieve a target stress. Rather, the target stress emerges from the various processes. We thus conclude that tensional homeostasis is not an end unto itself but a consequence of multiple cell and ECM processes in balance.

## Data accessibility

The datasets supporting this article have been uploaded as part of the electronic supplementary material.

## Authors' contributions

L.G. and V.H.B. contributed to the design of the study, writing of the MATLAB^{®} programs for the simulations, organization and analysis of the data. C.B.H. contributed to the planning of the study and coding. R.J.P. contributed to the data analysis. K.D.D. and Y.S. contributed to the planning of the study. All authors contributed to the writing of the manuscript.

## Competing interests

We have no competing interests.

## Funding

This work was supported by the National Science Foundation Graduate Research Fellowship (grant no. 00039202) to L.G. and by NSF grant no. CMMI 1300649 and NIH grant no. R01 EB005813-09S1.

## Acknowledgements

This work was carried out in part using computing resources at the University of Minnesota Supercomputing Institute and with recourses and the use of facilities at the Minneapolis VA Health Care System.

## Footnotes

One contribution of 19 to a theme issue ‘Integrated multiscale biomaterials experiment and modelling: towards function and pathology’.

- © 2015 The Author(s)