## Abstract

In soft tissues, large molecules such as proteoglycans trapped in the extracellular matrix (ECM) generate high levels of osmotic pressure to counter-balance external pressures. The semi-permeable matrix and fixed negative charges on these molecules serve to promote the swelling of tissues when there is an imbalance of molecular concentrations. Structural molecules, such as collagen fibres, form a network of stretch-resistant matrix, which prevents tissue from over-swelling and keeps tissue integrity. However, collagen makes little contribution to load bearing; the osmotic pressure in the ECM is the main contributor balancing external pressures. Although there have been a number of studies on tissue deformation, there is no rigorous analysis focusing on the contribution of the osmotic pressure in the ECM on the viscoelastic behaviour of soft tissues. Furthermore, most previous works were carried out based on the assumption of infinitesimal deformation, whereas tissue deformation is finite under physiological conditions. In the current study, a simplified mathematical model is proposed. Analytic solutions for solute distribution in the ECM and the free-moving boundary were derived by solving integro-differential equations under constant and dynamic loading conditions. Osmotic pressure in the ECM is found to contribute significantly to the viscoelastic characteristics of soft tissues during their deformation.

## 1. Introduction

Extracellular matrix (ECM) is ubiquitous in biological tissues and is composed of a mesh network of the fibrous protein collagen and, in some tissues, elastin bathed in an aqueous gel (Ayad *et al*. 1994). An individual collagen fibril is able to provide a substantial tensile strength when stretched, but makes little contribution in supporting compressive load (Winlove & Parker 1995). Structural investigation of collagenous meshwork reveals its regulating role in the transport of solute in connective tissues by excluding the volume accessible to solutes (Meyer *et al*. 1977; Bert *et al*. 1980). Owing to the exclusion volume, Maroudas *et al*. (1995) suggested a two-compartment model for the ECM with spaces accessible to macromolecules (i.e. inter-fibrillar space) and inaccessible to macromolecules (i.e. intra-fibrillar space). The major macromolecules in the ECM are proteoglycans (PGs) and other glycoproteins. They are trapped inside the ECM and are unable to permeate through the intra-fibrillar space. They are highly charged and attract a cloud of counter-ions around them. They generate osmotic pressure in the ECM and keep connective tissue ‘inflated’. The resultant osmotic pressure is restrained by the collagen network, which prevents the over-swelling of tissues and maintains tissue integrity (Comper & Laurent 1978). Osmotic pressure due to the PGs in the ECM is central to a wide range of physiological functions—enabling connective tissues such as articular cartilage to withstand high compressive loads and perform at very low frictional coefficient (Lai *et al*. 1991; Ateshian *et al*. 2003), influencing solute exchange through the interstitium (Maroudas 1979) and maintaining the structure and transparency of the cornea (Fatt 1978). A number of experimental and theoretical studies have been carried out to estimate the osmotic pressure by PGs in the ECM (Urban *et al*. 1979; Maroudas *et al*. 1991; Buschmann & Grodzinsky 1995; Ehrlich *et al*. 1998).

Many connective tissues are load-bearing, e.g. articular cartilage and intervertebral disc. Their deformation, mechanical response and metabolic regulation are affected by external mechanical loads. Different experimental protocols have been used to determine their intrinsic properties, such as the equilibrium modulus and the hydraulic permeability, which are employed in different theoretical models to simulate tissue deformation (Mow *et al*. 1984, 1999; Suh *et al*. 1995). In the biphasic model, the soft tissue is regarded as a mixture of a solid phase composed largely of ECM and a fluid phase consisting of water and its dissolved solutes (Mow *et al*. 1984). The frictional force between the two phases is used to explain the observed viscoelastic property of the tissue. In such a model, solid elastic stress exhibited by ECM and the pressurization of fluid phase work together to oppose the external load (Soltz & Ateshian 2000). The ratio of the contribution from these two factors varies with the loading condition.

Recognizing the importance of osmotic pressure arising from PGs and associated counter ions in the ECM, the ionic contribution was incorporated as a separate phase, leading to the triphasic theory (Lai *et al*. 1991) and quadraphasic theory (Huyghe & Janssen 1997). In these works, the role of large molecules, such as PGs, in the ECM is manifested by both the elastic stress of the ECM and the pressure of the fluid phase as well as the ionic contribution.

Most theoretical models have assumed infinitesimal tissue deformation in their analysis (Mow *et al*. 1984; Suh *et al*. 1995; Soltz & Ateshian 2000). Under physiological and experimental conditions, tissues usually undergo finite deformation with large local strain in the specimen observed in experiments (Holmes 1986). It is therefore important to incorporate finite deformation in a theoretical model. Furthermore, although the importance of the osmotic pressure of solutes in ECM on tissue deformation has been acknowledged qualitatively, there have been no rigorous theoretical studies that focus on the effect of the osmotic pressure generated by solutes in the ECM on external load bearing.

In the current study, a simplified mathematical model is developed to investigate effects of the osmotic pressure in ECM when the ECM is subjected to constant or dynamic loading conditions. To this end, mechanical stress of the matrix is neglected in order to explore the osmotic contribution to tissue deformation alone. Subject to the loading condition, the ECM can undergo finite deformation, which induces a moving boundary in the model. Effects of different external mechanical loads on the simplified system considering only the osmotic pressure are presented. Because of the idealized nature of the model, the results should not be considered to be predictive of physiological behaviour. They are only indicative of the complexity of the dynamic response of connective tissue when its osmotic pressure is considered. However, the theoretical results may be applicable to some simple *in vitro* experiments on the loading of connective tissue.

## 2. Theoretical model

A simplified one-dimensional model (as shown in figure 1) is considered in this study. The model is similar to the configuration in confined compression tests that are commonly adopted in experiments on connective tissue samples (Mow *et al*. 1984; Soltz & Ateshian 2000). In the model, the ECM containing large molecules (solutes) is confined to a container. All the walls of the container are impermeable to both the water and the solute, except the left-hand side of the container, which is separated from the bathing solution by a rigid, mass-less, semi-permeable membrane that allows water but not solute to pass. The membrane is free to move and mimics the porous filter used in confined compression tests. External mechanical loads are applied uniformly on the membrane. To separate the effect of the osmotic pressure from others, the meshwork is assumed to have zero-elasticity, i.e. the ECM acts as a connective tissue polysaccharide solution. Resistance from the system to external loads comes only from the osmotic pressure of the solute entrapped in the system. It should be noted that the highly charged macromolecules such as PGs and glycoproteins that make up the ground substances of most connective tissue are constrained within the tissue through entanglement with its fibril content and not by a semi-permeable membrane. Thus, the assumption of a massless semi-permeable membrane may not be such a bad model of real tissue.

For easy treatment of finite deformation in the model, the governing equations are set up in Lagrangian coordinates. Initial thickness of the model is *L* and the origin of the *x*-axis is set at the position of the membrane. At *t*=0, a uniform solute concentration *C*_{0} exists in the container. The membrane is pre-loaded and the initial load balances the osmotic pressure caused by *C*_{0}. When an external pressure is applied on the membrane, the position of the membrane changes with time and is denoted by *ξ*. At *t*=0, *ξ*=0. Mass conservation of the solute in the model gives:(2.1)where *C* is the concentration of the solute and *κ* is the solute diffusivity. In the current study, *κ* is assumed to be a constant, although it is generally a function of the concentration.

We assume that the membrane is so thin that its mass can be neglected and that the drag force caused by the water passing through it is also negligible. With these assumptions, it follows from local equilibrium that the osmotic pressure on the inner surface of the membrane, *π*(*t*), has to be equal to the applied stress on the outside of the membrane, *F*(*t*). Therefore, if the relationship between osmotic pressure and solute concentration is known, *π*=*Π*(*C*), the concentration of solute on the inner surface of the membrane is(2.2)where *Π*^{−1} is the inverse of the function *Π*. For the complex relationship between concentration and osmotic pressure typical of the highly charged solutes in connective tissue, it may not be possible to write this inverse function explicitly, but an implicit relationship can always be found due to the monotonically increasing nature of *Π*(*C*).

Owing to the unknown position of the membrane *ξ*, additional information is required to solve the problem. A Stefan condition is used in the vicinity of the moving membrane (Crank 1984)(2.3)where indicates the time derivative. It constitutes a second boundary condition at the moving membrane. On the other end of the container, the wall is impermeable,(2.4)The problem is completed with the initial condition:(2.5)

In the following analysis, *C*, *X*, *ξ* and *t* are non-dimensionalized using constant parameters, *C*_{0}, *L* and *κ*(2.6)where the tilde denotes non-dimensional variables. The non-dimensional equations are:(2.7)(2.8)(2.9)(2.10)(2.11)The tilde is omitted for writing convenience in the remaining part of this paper. Unless otherwise stated, all variables and parameters are non-dimensional.

## 3. Theoretical solutions

Equations (2.7)–(2.11) dictate the position of the moving boundary, i.e. the membrane and the distribution of the solute concentration in the solution. Moving boundaries occur in a class of problems, for example, multi-phase heat transfer (Kolodner 1956), oxygen diffusion (Hansen & Hougaard 1974) and diffusion-controlled growth and dissolution of precipitates (Enomoto & Atkinson 1993). Most theoretical solutions are restricted to the semi-infinite domain (e.g. Kolodner 1956; Crank 1984). Three major classes of numerical treatment have been used to tackle the moving-boundary problem. They are front-tracking, front-fixing and fixed domain (e.g. Crank 1984; Fazio 2001). In the front-tracking approach, the position of the moving boundary is computed explicitly by the numerical algorithm, whereas in the front-fixing or fixed domain approach the moving boundary is recovered *a posteriori* from the solution. These numerical methods often suffer from excessively long computation time and are unable to evaluate accurately the concentration gradient at the boundary (Elliot & Ockendon 1982; Enomoto & Atkinson 1993).

Enomoto & Atkinson (1993) developed a succinct method for diffusion-controlled growth of disordered interphase boundaries in a finite domain, based on theoretical results derived by Kolodner (1956) for the problem in a semi-infinite domain. The fixed boundary of the finite domain is treated by superimposing the solution from a fictitious boundary moving symmetrically with respect to the plane of the fixed boundary. The concentration variable along the container is analytically formulated as a function of the position of the moving membrane. The governing equations are thus reduced to the equation of the unknown position of the membrane *ξ*. To shorten expressions in the formulation and for easy comparison to the previous paper by Enomoto & Atkinson (1993), similar notation to that in Enomoto & Atkinson's work is used:(3.1)In the following analysis, equation (3.1) is used for the condition at the moving boundary instead of equation (2.10). This requires that:(3.2)From Kolodner (1956), the solution to equations (2.7)–(2.9), (2.11) and (3.1) over a semi-infinite domain (0<*X*<∞) can be expressed as:(3.3)where(3.4)(3.5)with *G* representing the Green function(3.6)

For the moving-boundary problem in a finite domain, which arises in the present model, the boundary condition on the fixed end of the container (i.e. equation (2.8)) is implemented by superimposing another solution for the mirror semi-infinite domain to equation (3.3). The detailed expression is given in §4 for different loading conditions on the membrane. The problem is now reduced to determine *ξ*, the position of the moving boundary, as a function of the time *t*. Following Kolodner's formulation, we consider a virtual domain on the opposite side of the moving boundary, the equation for the extended ‘whole’ domain can be derived as (for detailed derivations, see Kolodner 1956; Enomoto & Atkinson 1993):(3.7)This equation is an integro-differential equation, which determines the position of the moving boundary *ξ*. It has no analytical solution but is amenable to numerical treatment (Wei & Tao 1990). A similar method to that used by Enomoto & Atkinson (1993) is used to solve the equation, where the period of integration, *t*, is increased by a small amount iteratively. The moving velocity of the membrane is first obtained and is used to calculate the position of the membrane, *ξ*, in the next time-step. When *ξ* is solved, the concentration distribution can be easily calculated.

To initiate the iteration, time dependences of *ξ* and need to be provided. They are assumed to be the solutions in a semi-infinite domain when *t* is small. Enomoto & Atkinson (1993) give a detailed numerical treatment. When a constant load is applied, a fixed small initial time-step (*t*=0.001) is used which gives sufficiently accurate results. Under dynamic loading conditions, the loading frequency can be very high and results are sensitive to the choice of the initial time-step *t*. Therefore, for dynamic loads, the initial time-step *t* is selected as a hundredth of the loading period instead of a fixed small value. A simple linear relationship between the concentration of the solute and the osmotic pressure is used in following results (Tombs & Peacocke 1974). Since this relationship is only used to determine the instantaneous concentration at *ξ*, a more complicated, nonlinear relationship can be readily adopted in the model.

## 4. Results

Validation of the model has been carried out against results in Enomoto & Atkinson (1993) and good agreement has been confirmed. In the following results, different loading conditions are applied on the free-moving boundary. Mechanical responses of the system are present in terms of the solute concentration profiles and the change of the position of the moving boundary.

### (a) Constant load

A constant load indicates a step increase in the external pressure. For example, *f*(*t*)=1 indicates a step increase in the external pressure, which is equal to the osmotic pressure generated by solute with a concentration of *C*_{0}. When *Π*(*C*) is linear, as assumed here, it is convenient to define an intermediate constant parameter *Ω* in studies of constant loads of different magnitudes. *Ω*=1−*ξ* is the length of the solution in the container at equilibrium,(4.1)With equation (3.2),(4.2)From (3.3), concentration *C* can be expressed as:(4.3)Equations(3.7) and (4.3) lead to the integro-differential equation that determines the position of the moving boundary, *ξ*,(4.4)These equations are solved numerically using a finite difference method (Press *et al*. 1986).

Responses of the system to a constant load, *f*=1 (i.e. *Ω*=0.5) are presented in figure 2. In figure 2*a*, the velocity of the membrane is plotted against time (solid line). It is found that following the onset of the constant load, the membrane moves initially with a very high velocity because of the assumption that it is massless. This velocity decreases with time and approaches zero at *t*≈0.8. In comparison, for a semi-infinite domain, a theoretical solution for the velocity of the membrane exists and is shown as the dashed line in figure 2*a*. It is seen that the velocity of the membrane is initially the same as that in a finite domain. This is because of the diffusive nature of the process; the influence of the impermeable boundary at *X*=1 is insignificant until *t*∼0.2 when the flux of solute at the far boundary becomes significant. The effect of an impermeable wall on the membrane velocity becomes strong after *t*∼0.2, when the velocity of the membrane in a semi-finite domain decreases more quickly with time. For example, at *t*=0.8, the speed of the membrane in a semi-infinite domain is *ca* 0.5, whereas in a finite domain, the velocity has almost reached zero.

Concentration profiles of the solute inside the container at different times following the step increase in the pressure are presented in figure 2*b*. At *t*=0.01, a steep concentration gradient develops near the moving semi-permeable membrane, which is caused by the exudation of the water from the solution and the re-distribution of the solute in the container following the onset of the constant load. At this early stage, changes in the solute concentration are localized in an area close to the membrane. The elevated solute concentration in the vicinity of membrane generates an osmotic pressure, which counter-balances the increase in external pressure. As the solute diffuses down the concentration gradient away from the membrane into the inner regions of the container, the membrane moves further into the container as more water is exuded through the membrane. At *t*=0.1, concentration of the solute in the whole container has increased and the volume of the container is reduced by *ca* 27%. The concentration gradient inside the container falls gradually with the time as the system approaches its equilibrium state. At *t*=0.8, concentration of the solute in the container is almost uniform and the system approaches its equilibrium state. The volume of the container has been reduced by *ca* 50% by the application of this constant load.

### (b) Dynamic load

Dynamic loading of tissue specimens and cultured cells is commonly adopted in functional tissue engineering (e.g. Sah *et al*. 1992) and has been extensively modelled for the responses of the connective tissue to mechanical stimuli (e.g. Suh *et al*. 1995). In the current study, the external load is taken as a simple sinusoidal function,(4.5)where *ω* is the non-dimensional angular frequency of the loading pressure. *a* and *b* are constant parameters. In the present study, *a*=1 and *b*=0.5 are used as an example in the following results. From equation (3.2),(4.6)Accordingly, concentration of solute inside the container is given by:(4.7)Following the same procedure used for the constant load, the equation for the position of the membrane *ξ* is(4.8)The above equation is solved numerically using a finite difference method (Press *et al*. 1986).

The moving velocity of the membrane under cyclic loads with different loading frequencies is shown in figure 3. The biggest velocity occurs at the onset (*t*=0) of the external load. During the initial period of the time after the onset of the load, the velocity of the membrane decreases for all loading frequencies, although the magnitude of the external load is increasing. After the cyclic load reaches its peak value, the concentration of the solute near the membrane starts to decrease with the time, following the reduction in magnitude of the load. When the diffusive redistribution of the solute near the membrane cannot match the rate of reduction in the solute concentration there, water comes back into the container across the membrane. The volume of the container starts to expand and the direction of membrane movement alters. Steady-state oscillation is normally reached within 2–3 cycles of the load, where the velocity of the membrane alternates between the positive and the negative in a sinusoidal-like manner, following the cyclic load. The mechanical response of the system to cyclic loads is clearly different from that to a constant load. During the initial period of time, the magnitude of the membrane velocity under cyclic loads is higher than that under a constant load. This is understandable since the magnitude of the load that we are applying is higher during this period for the cyclic load. Another major difference lies in the direction of the movement of the membrane. Under constant load, the membrane moves in one direction only towards the end of container. Under cyclic load, the velocity of the membrane not only changes in its magnitude but also alternates the directions. Higher loading frequency results in bigger magnitude of the velocity of the membrane (see figure 6).

Distribution of the solute concentration in the container under cyclic load exhibits different profiles from those under constant load. The concentration is no longer always bigger near the membrane as observed in the container under constant load. In figure 4, solute concentration distributions in the container at different times are presented when the system has reached steady oscillation state. As the magnitude of the load increases, the membrane moves into the container and water is exuded. Higher solute concentration exists near the membrane and solute diffuses into the inner region of the container. When the external load starts to decrease, the corresponding concentration of the solute near the membrane is reduced instantaneously and water is imbibed into the container. This reduction could lower the solute concentration in the vicinity of the membrane to less than those in the inner regions of the container.

Figure 5 shows the temporal changes of the velocity (figure 5*a*) and the displacement (figure 5*b*) of the membrane in steady oscillation state. The time axis is converted to the number of loading cycles, i.e. *tω*/2*π*. It is found that the loading frequency plays an important role in the mechanical response of the system to cyclic loads. At a low loading frequency (e.g. *ω*=1/8), the membrane moves at a very slow velocity (e.g. peak value of the velocity is *ca* 0.015). There is not much movement across the membrane. The slow variation in the loading magnitude with the time is largely met by the redistribution (diffusion) of the solute inside the container. If we compare the phase difference between the velocity of the membrane and the external load, we find the peak value of the membrane velocity is *ca* 0.656*π* ahead of the external load. Although the velocity of the membrane is small, the low frequency of the load provides the membrane with enough time for displacement. It is seen in figure 5*b* that the membrane moves between *x*=0.41 and *x*=0.65 at *ω*=1/8. As the loading frequency is increased to *ω*=1, the peak velocity of the membrane increases to 0.28. The phase difference between the velocity and the external load is *ca* 0.561*π*. Movement of the membrane is between *x*=0.37 and *x*=0.61. At even higher loading frequency, *ω*=8, the peak velocity of the membrane increases to 0.82. The phase difference between the membrane velocity and the external load is only *ca* 0.339*π*. The membrane moves between *x*=0.39 and *x*=0.59. The range of the movement of the membrane at *ω*=8 (i.e. 0.20) is smaller than those at *ω*=1 or 1/8 (i.e. 0.24). Increases in the loading frequency lead to higher membrane velocity, smaller range of the displacement of the membrane and reduction in the phase difference between velocity of the membrane and the load.

Effects of the loading frequency on the peak velocity of the membrane and the range of the movement of the membrane in steady oscillation state are presented in figure 6. It is found that the peak velocity of the membrane increases rapidly (and almost proportionally) with the loading frequency. On the other hand, the range of the movement of the membrane is rather insensitive to changes in the loading frequency. It reaches maximum at *ω*≈1 and then decreases with the loading frequency. It is noteworthy that the membrane velocity curves in figure 5 are not strictly sinusoidal. The curves skew strongly at low loading frequencies, e.g. *ω*=1. The curves of membrane position with the time in figure 5 are more sinusoidal-like, in comparison, over the range of frequencies investigated. These are in line with the observed behaviour of soft tissue samples under cyclic loads (Suh *et al*. 1995).

## 5. Discussion

This study uses a simple theoretical model to explore the effect of osmotic pressure on the dynamics of the ECM under constant and cyclic external loads. All other factors that could have contributed to load bearing in the ECM are neglected. This serves to highlight the effect of the osmotic pressure in the ECM, but also brings clear limitations to the direct comparison of the model to experimental results. The membrane is assumed to be rigid and massless. There is no drag between the water and the membrane when water moves through the membrane. These assumptions are necessary to relate external load to the osmotic pressure in the vicinity of the membrane. In real situations, drag force between the fluid and the membrane may contribute to the force balance on the membrane and alter its movement.

In this study, the membrane has an initial load that counter-balances the osmotic pressure generated by the solute concentration *C*_{0} in the container. When additional load is applied on the membrane, the system changes (e.g. movement of the membrane and re-distribution of the solute in the container) in response to the additional load applied. This has similarity to experiments, where zero-tension of the collage meshwork is first achieved before the deformation of the connective tissue is investigated. For example, when a cartilage specimen is soaked in a hypertonic solution, the collagen mesh could ‘crimp’ and be close to zero-stress when the osmotic pressure of PGs in the cartilage is balanced by the osmotic pressure of the bathing solution. The dimension of cartilage samples has been observed to decrease with increasing NaCl concentration in experiments (Lai *et al*. 1991).

Osmotic pressure in the ECM is generated mainly from the negative fixed charges on glycosaminoglycan chains of PGs (Maroudas 1968; Comper & Laurent 1978). In the current model, osmotic pressure is directly related to the concentration of the PGs, since counter ions are closely associated with the fixed charges on PGs and their concentration is related to the fixed charge density (FCD) on PGs. The theoretical model can be used to model osmotic pressure by neutral molecules, as well as charged molecules. With charged molecules the FCD can be chosen as the variable instead of the concentration of the molecule (Urban *et al*. 1979).

The relationship between the osmotic pressure and the solute concentration for macromolecules, *Π*(*C*), is often nonlinear (Maroudas *et al*. 1995; Ehrlich *et al*. 1998). The model can be readily adapted to accommodate this complexity since *Π*(*C*) is only used to determine the osmotic pressure at the semi-permeable membrane. Furthermore, it can be used to explore the relationship between the osmotic pressure and solute concentration.

## 6. Conclusion

A simplified mathematical model has been developed to examine the contribution of the osmotic pressure generated by solutes in the ECM to external load bearing. The model should not be seen as a realistic description of physiological processes; they will involve more complex, three-dimensional geometries and contributions from other load-bearing components of the tissue. Instead, it should be seen as an indication of the important effect that the osmotic pressure, frequently overlooked, can have on connective tissue mechanics. The simplified model has been chosen to highlight these effects.

Although it is simple, the model still entails complex mathematical analysis. Theoretical results on the deformation of the tissue and concentration profiles of the solute in extracellular fluid have been achieved by solving integro-differentiation equations. Under constant loads, the membrane yields towards a new equilibrium state and the movement of the membrane exhibits behaviour typical of viscoelastic processes although it is generated by diffusion. The concentration of the solute in the extracellular fluid increases with the time and reaches a new elevated uniform distribution at equilibrium. Under cyclic loads, movement of the membrane reaches a steady oscillation state after the load is applied for a few cycles. Higher loading frequency results in bigger velocity of the membrane, but has little effect on the displacement of the membrane. The simplified model neglects other contributions to load bearing in ECM but provides insight into the role of osmotic pressure generated by macromolecule in the ECM under different loading conditions. Although the direct application of the model to physiological processes is problematic, it could be used to interpret some *in vitro* experiments, where one-dimensional, simple loading is applied to tissue constrained within an impermeable container.

## Editors' note

Please see also related communications in this focussed issue by McDougall *et al.* (2006) and Schaller & Meyer-Hermann (2006).

## Acknowledgments

The authors would like to thank Professor Colin Atkinson, FRS, Imperial College, London, for helpful discussions.

## Footnotes

One contribution of 15 to a Theme Issue ‘Biomathematical modelling II’.

- © 2006 The Royal Society