Sensitivity of cartilage mechanical behaviour to spatial variations in material properties.

Articular cartilage tissue exhibits a spatial dependence in material properties that govern mechanical behaviour. A mathematical model of cartilage tissue under one dimensional confined compression testing is developed for normal tissue that takes account of these variations in material properties. Modifications to the model representative of a selection of mechanisms driving osteoarthritic cartilage are proposed, allowing application of the model to both physiological and pathophysiological, osteoarthritic tissue. Incorporating spatial variations into the model requires the specification of more parameters than are required in the absence of these variations. A global sensitivity analysis of these parameters is implemented to identify the dominant mechanisms of mechanical response, in normal and osteoarthritic cartilage tissue, to both static and dynamic loading. The most sensitive parameters differ between dynamic and static mechanics of the cartilage, and also differ between physiological and osteoarthritic pathophysiological cartilage. As a consequence changes in cartilage mechanics in response to alterations in cartilage structure are predicted to be contingent on the nature of loading and the health, or otherwise, of the cartilage. In particular the mechanical response of cartilage, especially deformation, is predicted to be much more sensitive to cartilage stiffness in the superficial zone given the onset of osteoarthritic changes to material properties, such as superficial zone increases in permeability and reductions in fixed charge. In turn this indicates that any degenerative changes in the stiffness associated with the superficial cartilage collagen mesh are amplified if other elements of osteoarthritic disease are present, which provides a suggested mechanism-based explanation for observations that the range of mechanical parameters representative of normal and osteoarthritic tissue can overlap substantially.

Modelling articular cartilage: the relative motion of two adjacent poroviscoelastic layers.

In skeletal joints two layers of adjacent cartilage are often in relative motion. The individual cartilage layers are often modelled as a poroviscoelastic material. To model the relative motion, noting the separation of scales between the pore level and the macroscale, a homogenization based on multiple scale asymptotic analysis has been used in this study to derive a macroscale model for the relative translation of two poroviscoelastic layers separated by a very thin layer of fluid. In particular the fluid layer thickness is essentially zero at the macroscale so that the two poroviscoelastic layers are effectively in contact and their interaction is captured in the derived model via a set of interfacial conditions, including a generalization of the Beavers-Joseph condition at the interface between a viscous fluid and a porous medium. In the simplifying context of a uniform geometry, constant fixed charge density, a Newtonian interstitial fluid and a viscoelastic scaffold, modelled via finite deformation theory, we present preliminary simulations that may be used to highlight predictions for how oscillatory relative movement of cartilage under load influences the peak force the cartilage experiences and the extent of the associated deformations. In addition to highlighting such cartilage mechanics, the systematic derivation of the macroscale models will enable the study of how nanoscale cartilage physics, such as the swelling pressure induced by fixed charges, manifests in cartilage mechanics at much higher lengthscales.

Model reduction for initial value ODEs.

Many physical phenomena in biology and physiology are described by mathematical models that comprise a system of initial value ordinary differential equations. Each differential equation may often be written as the sum of several terms, where each term represents a different physical entity. A wide range of techniques, ranging from heuristic observation to mathematically rigorous asymptotic analysis, may be used to simplify these equations allowing the identification of the key phenomena responsible for a given observed behaviour. In this study we extend an algorithm for automatically simplifying systems of initial value ordinary differential equations (Whiteley (2010)) that is based on a posteriori analysis of the full system of equations. Our extensions to the algorithm make the following contributions: (i) each equation in a system of differential equations may be written as a finite sum of contributions (including the derivative term), and any one of these terms may be neglected (if it is appropriate to do so) in the simplified model; and (ii) a simplified model is generated that allows accurate prediction of one or more components of the solution at all times. These extensions are illustrated using examples drawn from enzyme kinetics and cardiac electrophysiology.

Modelling the inclusion of swelling pressure in a tissue level poroviscoelastic model of cartilage deformation.

Swelling pressure in the interstitial fluid within the pores of cartilage tissue is known to have a significant effect on the rheology of cartilage tissue. The swelling pressure varies rapidly within thin regions inside pores known as Debye layers, caused by the presence of fixed charge, as observed in cartilage. Tissue level calculation of cartilage deformation therefore requires resolution of three distinct spatial scales: the Debye lengthscale within individual pores; the lengthscale of an individual pore; and the tissue lengthscale. We use asymptotics to construct a leading order approximation to the swelling pressure within pores, allowing the swelling pressure to be systematically included within a fluid-solid interaction model at the level of pores in cartilage. We then use homogenization to derive tissue level equations for cartilage deformation that are very similar to those governing the finite deformation of a poroviscoelastic body. The equations derived permit the spatial variations in porosity and electric charge that occur in cartilage tissue. Example solutions are then used to confirm the plausibility of the model derived and to consider the impact of fixed charge heterogeneity, illustrating that local fixed charge loss is predicted to increase deformation gradients under confined compression away from, rather than at, the site of loss.

An evaluation of some assumptions underpinning the bidomain equations of electrophysiology.

Tissue level cardiac electrophysiology is usually modelled by the bidomain equations, or the monodomain simplification of the bidomain equations. One assumption made when deriving the bidomain equations is that both the intracellular and extracellular spaces are in electrical equilibrium. This assumption neglects the disturbance of this equilibrium in thin regions close to the cell membrane known as Debye layers. We first demonstrate that the governing equations at the cell, or microscale, level may be adapted to take account of these Debye layers with little additional complexity, provided the permittivity within the Debye layers satisfies certain conditions that are believed to be satisfied for biological cells. We then homogenize the microscale equations using a technique developed for an almost periodic microstructure. Cardiac tissue is usually modelled as sheets of cardiac fibres stacked on top of one another. A common assumption is that an orthogonal coordinate system can be defined at each point of cardiac tissue, where the first axis is in the fibre direction, the second axis is orthogonal to the first axis but lies in the sheet of cardiac fibres and the third axis is orthogonal to the cardiac sheet. It is assumed further that both the intracellular and extracellular conductivity tensors are diagonal with respect to these axes and that the diagonal entries of these tensors are constant across the whole tissue. Using the homogenization technique we find that this assumption is usually valid for cardiac tissue, but highlight situations where the assumption may not be valid.

Identifying chondrogenesis strategies for tissue engineering of articular cartilage.

A key step in the tissue engineering of articular cartilage is the chondrogenic differentiation of mesenchymal stem cells (MSCs) into chondrocytes (native cartilage cells). Chondrogenesis is regulated by transforming growth factor-ß (TGF-ß), a short-lived cytokine whose effect is prolonged by storage in the extracellular matrix. Tissue engineering applications aim to maximise the yield of differentiated MSCs. Recent experiments involve seeding a hydrogel construct with a layer of MSCs lying below a layer of chondrocytes, stimulating the seeded cells in the construct from above with exogenous TGF-ß and then culturing it in vitro. To investigate the efficacy of this strategy, we develop a mathematical model to describe the interactions between MSCs, chondrocytes and TGF-ß. Using this model, we investigate the effect of varying the initial concentration of TGF-ß, the initial densities of the MSCs and chondrocytes, and the relative depths of the two layers on the long-time composition of the tissue construct.

The combined impact of tissue heterogeneity and fixed charge for models of cartilage: the one-dimensional biphasic swelling model revisited.

Articular cartilage is a complex, anisotropic, stratified tissue with remarkable resilience and mechanical properties. It has been subject to extensive modelling as a multiphase medium, with many recent studies examining the impact of increasing detail in the representation of this tissue's fine scale structure. However, further investigation of simple models with minimal constitutive relations can nonetheless inform our understanding at the foundations of soft tissue simulation. Here, we focus on the impact of heterogeneity with regard to the volume fractions of solid and fluid within the cartilage. Once swelling pressure due to cartilage fixed charge is also present, we demonstrate that the multiphase modelling framework is substantially more complicated, and thus investigate this complexity, especially in the simple setting of a confined compression experiment. Our findings highlight the importance of locally, and thus heterogeneously, approaching pore compaction for load bearing in cartilage models, while emphasising that such effects can be represented by simple constitutive relations. In addition, simulation predictions are observed for the sensitivity of stress and displacement in the cartilage to variations in the initial state of the cartilage and thus the details of experimental protocol, once the tissue is heterogeneous. These findings are for the simplest models given only heterogeneity in volume fractions and swelling pressure, further emphasising that the complex behaviours associated with the interaction of volume fraction heterogeneity and swelling pressure are likely to persist for simulations of cartilage representations with more fine-grained structural detail of the tissue.

Representation of Multiple Cellular Phenotypes Within Tissue-Level Simulations of Cardiac Electrophysiology.

Distinct electrophysiological phenotypes are exhibited by biological cells that have differentiated into particular cell types. The usual approach when simulating the cardiac electrophysiology of tissue that includes different cell types is to model the different cell types as occupying spatially distinct yet coupled regions. Instead, we model the electrophysiology of well-mixed cells by using homogenisation to derive an extension to the commonly used monodomain or bidomain equations. These new equations permit spatial variations in the distribution of the different subtypes of cells and will reduce the computational demands of solving the governing equations. We validate the homogenisation computationally, and then use the new model to explain some experimental observations from stem cell-derived cardiomyocyte monolayers.

Predictive Mathematical Models for the Spread and Treatment of Hyperoxia-induced Photoreceptor Degeneration in Retinitis Pigmentosa.

Purpose: To determine whether the oxygen toxicity hypothesis can explain the distinctive spatio-temporal patterns of retinal degeneration associated with human retinitis pigmentosa (RP) and to predict the effects of antioxidant and trophic factor treatments under this hypothesis. Methods: Three mathematical models were derived to describe the evolution of the retinal oxygen concentration and photoreceptor density over time. The first model considers only hyperoxia-induced degeneration, while the second and third models include mutation-induced rod and cone loss respectively. The models were formulated as systems of partial differential equations, defined on a two-dimensional domain spanning the region between the foveal center and the ora serrata, and were solved numerically using the finite element method. Results: The mathematical models recapitulate patterns of retinal degeneration which involve preferential loss of photoreceptors in the parafoveal/perifoveal and far-peripheral retina, while those which involve a preferential loss of midperipheral photoreceptors cannot be reproduced. Treatment with antioxidants or trophic factors is predicted to delay, halt, or partially reverse retinal degeneration, depending upon the strength and timing of treatment and disease severity. Conclusions: The model simulations indicate that while the oxygen toxicity hypothesis is sufficient to explain some of the patterns of retinal degeneration observed in human RP, additional mechanisms are necessary to explain the full range of behaviors. The models further suggest that antioxidant and trophic factor treatments have the potential to reduce hyperoxia-induced disease severity and that, where possible, these treatments should be targeted at retinal regions with low photoreceptor density to maximize their efficacy.

Potential for noninvasive assessment of lung inhomogeneity using highly precise, highly time-resolved measurements of gas exchange.

Inhomogeneity in the lung impairs gas exchange and can be an early marker of lung disease. We hypothesized that highly precise measurements of gas exchange contain sufficient information to quantify many aspects of the inhomogeneity noninvasively. Our aim was to explore whether one parameterization of lung inhomogeneity could both fit such data and provide reliable parameter estimates. A mathematical model of gas exchange in an inhomogeneous lung was developed, containing inhomogeneity parameters for compliance, vascular conductance, and dead space, all relative to lung volume. Inputs were respiratory flow, cardiac output, and the inspiratory and pulmonary arterial gas compositions. Outputs were expiratory and pulmonary venous gas compositions. All values were specified every 10 ms. Some parameters were set to physiologically plausible values. To estimate the remaining unknown parameters and inputs, the model was embedded within a nonlinear estimation routine to minimize the deviations between model and data for CO2, O2, and N2 flows during expiration. Three groups, each of six individuals, were studied: young (20-30 yr); old (70-80 yr); and patients with mild to moderate chronic obstructive pulmonary disease (COPD). Each participant undertook a 15-min measurement protocol six times. For all parameters reflecting inhomogeneity, highly significant differences were found between the three participant groups ( P < 0.001, ANOVA). Intraclass correlation coefficients were 0.96, 0.99, and 0.94 for the parameters reflecting inhomogeneity in deadspace, compliance, and vascular conductance, respectively. We conclude that, for the particular participants selected, highly repeatable estimates for parameters reflecting inhomogeneity could be obtained from noninvasive measurements of respiratory gas exchange. NEW & NOTEWORTHY This study describes a new method, based on highly precise measures of gas exchange, that quantifies three distributions that are intrinsic to the lung. These distributions represent three fundamentally different types of inhomogeneity that together give rise to ventilation-perfusion mismatch and result in impaired gas exchange. The measurement technique has potentially broad clinical applicability because it is simple for both patient and operator, it does not involve ionizing radiation, and it is completely noninvasive.

Inducing chondrogenesis in MSC/chondrocyte co-cultures using exogenous TGF-ß: a mathematical model.

The differentiation of mesenchymal stem cells (MSCs) into chondrocytes (native cartilage cells), or chondrogenesis, is a key step in the tissue engineering of articular cartilage, where the motility and high proliferation rate of MSCs used as seed cells are exploited. Chondrogenesis is regulated by transforming growth factor-beta (TGF-ß), a short-lived cytokine whose effect is prolonged by storage in the extracellular matrix. Tissue engineering applications require the complete differentiation of an initial population of MSCs, and two common strategies used to achieve this in vitro are (1) co-culture the MSCs with chondrocytes, which constitutively produce TGF-ß; or (2) add exogenous TGF-ß. To investigate these strategies we develop an ordinary differential equation model of the interactions between TGF-ß, MSCs and chondrocyte. Here the dynamics of TGF-ß are much faster than those of the cell processes; this difference in time-scales is exploited to simplify subsequent model analysis. Using our model we demonstrate that under strategy 1 complete chondrogenesis will be induced if the initial proportion of chondrocytes exceeds a critical value. Similarly, under strategy 2 we find that there is a critical concentration of exogenous TGF-ß above which all MSCs will ultimately differentiate. Finally, we use the model to demonstrate the potential advantages of adopting a hybrid strategy where exogenous TGF-ß is added to a co-culture of MSCs and chondrocytes, as compared to using either strategy 1 or 2 in isolation.

Mathematical modelling of cell layer growth in a hollow fibre bioreactor.

Generating autologous tissue grafts of a clinically useful volume requires efficient and controlled expansion of cell populations harvested from patients. Hollow fibre bioreactors show promise as cell expansion devices, owing to their potential for scale-up. However, further research is required to establish how to specify appropriate hollow fibre bioreactor operating conditions for expanding different cell types. In this study we develop a simple model for the growth of a cell layer seeded on the outer surface of a single fibre in a perfused hollow fibre bioreactor. Nutrient-rich culture medium is pumped through the fibre lumen and leaves the bioreactor via the lumen outlet or passes through the porous fibre walls and cell layer, and out via ports on the outer wall of the extra-capillary space. Stokes and Darcy equations for fluid flow in the fibre lumen, fibre wall, cell layer and extra-capillary space are coupled to reaction-advection-diffusion equations for oxygen and lactate transport through the bioreactor, and to a simple growth law for the evolution of the free boundary of the cell layer. Cells at the free boundary are assumed to proliferate at a rate that increases with the local oxygen concentration, and to die and detach from the layer if the local fluid shear stress or lactate concentration exceed critical thresholds. We use the model to predict operating conditions that maximise the cell layer growth for different cell types. In particular, we predict the optimal flow rate of culture medium into the fibre lumen and fluid pressure imposed at the lumen outlet for cell types with different oxygen demands and fluid shear stress tolerances, and compare the growth of the cell layer when the exit ports on the outside of the bioreactor are open with that when they are closed. Model simulations reveal that increasing the inlet flow rate and outlet fluid pressure increases oxygen delivery to the cell layer and, therefore, the growth rate of cells that are tolerant to high shear stresses, but may be detrimental for shear-sensitive cells. The cell layer growth rate is predicted to increase, and be less sensitive to the lactate tolerance of the cells, when the exit ports are opened, as the radial flow through the bioreactor is enhanced and the lactate produced by the cells cleared more rapidly from the cell layer.

Evaluation of the growth environment of a hydrostatic force bioreactor for preconditioning of tissue-engineered constructs.

Bioreactors have been widely acknowledged as valuable tools to provide a growth environment for engineering tissues and to investigate the effect of physical forces on cells and cell-scaffold constructs. However, evaluation of the bioreactor environment during culture is critical to defining outcomes. In this study, the performance of a hydrostatic force bioreactor was examined by experimental measurements of changes in dissolved oxygen (O2), carbon dioxide (CO2), and pH after mechanical stimulation and the determination of physical forces (pressure and stress) in the bioreactor through mathematical modeling and numerical simulation. To determine the effect of hydrostatic pressure on bone formation, chick femur skeletal cell-seeded hydrogels were subjected to cyclic hydrostatic pressure at 0-270 kPa and 1 Hz for 1 h daily (5 days per week) over a period of 14 days. At the start of mechanical stimulation, dissolved O2 and CO2 in the medium increased and the pH of the medium decreased, but remained within human physiological ranges. Changes in physiological parameters (O2, CO2, and pH) were reversible when medium samples were placed in a standard cell culture incubator. In addition, computational modeling showed that the distribution and magnitude of physical forces depends on the shape and position of the cell-hydrogel constructs in the tissue culture format. Finally, hydrostatic pressure was seen to enhance mineralization of chick femur skeletal cell-seeded hydrogels.

Optimising cell aggregate expansion in a perfused hollow fibre bioreactor via mathematical modelling.

The need for efficient and controlled expansion of cell populations is paramount in tissue engineering. Hollow fibre bioreactors (HFBs) have the potential to meet this need, but only with improved understanding of how operating conditions and cell seeding strategy affect cell proliferation in the bioreactor. This study is designed to assess the effects of two key operating parameters (the flow rate of culture medium into the fibre lumen and the fluid pressure imposed at the lumen outlet), together with the cell seeding distribution, on cell population growth in a single-fibre HFB. This is achieved using mathematical modelling and numerical methods to simulate the growth of cell aggregates along the outer surface of the fibre in response to the local oxygen concentration and fluid shear stress. The oxygen delivery to the cell aggregates and the fluid shear stress increase as the flow rate and pressure imposed at the lumen outlet are increased. Although the increased oxygen delivery promotes growth, the higher fluid shear stress can lead to cell death. For a given cell type and initial aggregate distribution, the operating parameters that give the most rapid overall growth can be identified from simulations. For example, when aggregates of rat cardiomyocytes that can tolerate shear stresses of up to 0:05 Pa are evenly distributed along the fibre, the inlet flow rate and outlet pressure that maximise the overall growth rate are predicted to be in the ranges 2.75 x 10(-5) m(2) s(-1) to 3 x 10(-5) m(2) s(-1) (equivalent to 2.07 ml min(-1) to 2.26 ml min(-1)) and 1.077 x 10(5) Pa to 1.083 x 10(5) Pa (or 15.6 psi to 15.7 psi) respectively. The combined effects of the seeding distribution and flow on the growth are also investigated and the optimal conditions for growth found to depend on the shear tolerance and oxygen demands of the cells.

Automatic simplification of systems of reaction-diffusion equations by a posteriori analysis.

Many mathematical models in biology and physiology are represented by systems of nonlinear differential equations. In recent years these models have become increasingly complex in order to explain the enormous volume of data now available. A key role of modellers is to determine which components of the model have the greatest effect on a given observed behaviour. An approach for automatically fulfilling this role, based on a posteriori analysis, has recently been developed for nonlinear initial value ordinary differential equations [J.P. Whiteley, Model reduction using a posteriori analysis, Math. Biosci. 225 (2010) 44-52]. In this paper we extend this model reduction technique for application to both steady-state and time-dependent nonlinear reaction-diffusion systems. Exemplar problems drawn from biology are used to demonstrate the applicability of the technique.

Modelling the effect of gap junctions on tissue-level cardiac electrophysiology.

When modelling tissue-level cardiac electrophysiology, a continuum approximation to the discrete cell-level equations, known as the bidomain equations, is often used to maintain computational tractability. Analysing the derivation of the bidomain equations allows us to investigate how microstructure, in particular gap junctions that electrically connect cells, affect tissue-level conductivity properties. Using a one-dimensional cable model, we derive a modified form of the bidomain equations that take gap junctions into account, and compare results of simulations using both the discrete and continuum models, finding that the underlying conduction velocity of the action potential ceases to match up between models when gap junctions are introduced at physiologically realistic coupling levels. We show that this effect is magnified by: (i) modelling gap junctions with reduced conductivity; (ii) increasing the conductance of the fast sodium channel; and (iii) an increase in myocyte length. From this, we conclude that the conduction velocity arising from the bidomain equations may not be an accurate representation of the underlying discrete system. In particular, the bidomain equations are less likely to be valid when modelling certain diseased states whose symptoms include a reduction in gap junction coupling or an increase in myocyte length.

A numerical guide to the solution of the bi-domain equations of cardiac electrophysiology.

Simulation of cardiac electrical activity using the bi-domain equations can be a massively computationally demanding problem. This study provides a comprehensive guide to numerical bi-domain modelling. Each component of bi-domain simulations--discretization, ODE-solution, linear system solution, and parallelization--is discussed, and previously-used methods are reviewed, new methods are proposed, and issues which cause particular difficulty are highlighted. Particular attention is paid to the choice of stimulus currents, compatibility conditions for the equations, the solution of singular linear systems, and convergence of the numerical scheme.

Modeling alveolar volume changes during periodic breathing in heterogeneously ventilated lungs.

A simplified model of periodic breathing, proposed by Whiteley et al. (Math. Med. Biol. 20:205-224, 2003), describes a non-uniform breathing pattern for a lung with an inhomogeneous gas distribution, such as that observed in some subjects suffering from respiratory disease. This model assumes a constant alveolar volume, and predicts incidence of irregular breathing caused by small, poorly ventilated regions of the lung. Presented here is an extension to this work which, by allowing variable lung volume, facilitates the investigation of pulmonary collapse in poorly ventilated compartments. A weakness of the original model is that a very small alveolar volume is required for periodic breathing to occur. The model presented within, which removes the assumption of constant compartment volume and allows alveolar volume to vary with time, predicts periodic breathing at higher, more realistic alveolar volumes. Furthermore, the predicted oscillations in ventilation match experimental data more closely. Thus the model that allows for alveolar collapse has improved upon these earlier results, and establishes a theoretical link between periodic breathing and atelectasis.

Model reduction using a posteriori analysis.

Mathematical models in biology and physiology are often represented by large systems of non-linear ordinary differential equations. In many cases, an observed behaviour may be written as a linear functional of the solution of this system of equations. A technique is presented in this study for automatically identifying key terms in the system of equations that are responsible for a given linear functional of the solution. This technique is underpinned by ideas drawn from a posteriori error analysis. This concept has been used in finite element analysis to identify regions of the computational domain and components of the solution where a fine computational mesh should be used to ensure accuracy of the numerical solution. We use this concept to identify regions of the computational domain and components of the solution where accurate representation of the mathematical model is required for accuracy of the functional of interest. The technique presented is demonstrated by application to a model problem, and then to automatically deduce known results from a cell-level cardiac electrophysiology model.

Solving the coupled system improves computational efficiency of the bidomain equations.

The bidomain equations are frequently used to model the propagation of cardiac action potentials across cardiac tissue. At the whole organ level, the size of the computational mesh required makes their solution a significant computational challenge. As the accuracy of the numerical solution cannot be compromised, efficiency of the solution technique is important to ensure that the results of the simulation can be obtained in a reasonable time while still encapsulating the complexities of the system. In an attempt to increase efficiency of the solver, the bidomain equations are often decoupled into one parabolic equation that is computationally very cheap to solve and an elliptic equation that is much more expensive to solve. In this study, the performance of this uncoupled solution method is compared with an alternative strategy in which the bidomain equations are solved as a coupled system. This seems counterintuitive as the alternative method requires the solution of a much larger linear system at each time step. However, in tests on two 3-D rabbit ventricle benchmarks, it is shown that the coupled method is up to 80% faster than the conventional uncoupled method-and that parallel performance is better for the larger coupled problem.