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.


Introduction
Articular cartilage is a thin layer of dense tissue found on the ends of bones in synovial joints. It prevents contact between the ends of the bones and acts as a lubricating surface between them [1,2].
Cartilage is composed of a combination of native cartilage cells (or chondrocytes) and extracellular matrix  : Schematic cartoon of the zonated, depth-dependent structure of articular cartilage. In a synovial joint the cartilage lies between the synovial fluid and the subchondral bone, and is typically described as being divided into three zones. The highest density of both collagen and chondrocytes is in the superficial zone, and here both are aligned in the direction of the surface of the cartilage. In the middle zone the cells are more rounded and the alignment of the collagen is less organised, with the density of both lower than in the superficial zone. In the deep zone the cells are organised in column-like structures, with both these columns and the collagen aligned perpendicular to the subchondral bone [2].
components, principally collagen and proteoglycans. The zonated structure and arrangement of both the cells and the matrix components endow cartilage with properties that make it particularly well suited to withstanding mechanical stresses in load bearing joints, such as the knee (a schematic cross-section of this zonated structure is shown in Figure 1). Articular cartilage has a low capacity to regenerate due to its avascular structure and the low motility and proliferation rates of the chondrocytes. Thus, any degeneration of cartilage is chronic and can lead to severe conditions, such as osteoarthritis, which typically require some form of surgical intervention. A proposed alternative to current joint replacement approaches involves implanting artificially engineered cartilage to replace the damaged tissue.
Ideally such an implant should mimic the biomechanical function of the natural tissue and have a similar distribution of cells and matrix components. One promising effort to produce such an implant involves seeding cells in a hydrogel construct reinforced with a lattice of 3D-printed polymer fibres, which is then cultured in a bioreactor under chemical and mechanical stimuli [3]. The seeded cells may be either chondrocytes harvested from natural cartilage, or mesenchymal stem cells (MSCs), which differentiate into chondrocytes in response to externally imposed stimuli, or a combination of both cell types. The relevant matrix components are synthesised and maintained by chondrocytes, so a zonated distribution of these components requires a zonated distribution of the chondrocytes. As chondrocytes have very low motility this distribution may be achieved either by seeding chondrocytes in a zonated fashion, or by exploiting the high motility of the MSCs and stimulating them to move throughout the construct before they differentiate.
A further consideration is the low proliferation rate of chondrocytes which means that many cells must be harvested from other sources to seed the construct; such harvesting is practically difficult to achieve.
By contrast, MSCs are highly proliferative, so fewer cells need be harvested; the challenges in their use are to bias their differentiation into chondrocytes (chondrogenesis) and then maintain a stable phenotype [4] (we note that MSCs may also differentiate into bone or muscle cells, for instance). Chemical stimulation can be provided to the seeded cells, for instance in the form of fetal bovine serum (FBS), platelet-derived growth factor (PDGF), or transforming growth factor-β (TGF-β), to promote proliferation, differentiation and/or chemotactic movement of the MSCs [5,6,7]. Chemical stimulation can be combined with mechanical stimulation, for instance by using a piston to apply a load to the hydrogel construct [8], to promote cell differentiation via mechanotransduction.
Identifying appropriate conditions under which to culture artificial constructs to produce a specific endproduct is technically challenging. The cells are typically highly sensitive to small changes in the environment, and therefore a small change in the concentration of a particular growth factor may have a large effect.
Similarly, for mechanical loading a balance needs to be struck between applying enough loading to induce mechanotransductive effects, but not too much that the structure of the hydrogel is compromised. Additionally, in the context of tissue engineering differentiated cells may dedifferentiate, and so may require regular stimulation to prevent their loss of phenotype. Producing a desired outcome involves carefully balancing the application of the different stimuli and the responses they induce in the cells. The modelling work of this study aims to address this challenging question. We develop a mathematical model to investigate how to drive a population of MSCs to chondrogenesis with a focus on the response of MSCs to a single growth factor in particular, transforming growth factor β (TGF-β).
TGF-β is particularly important in a number of contexts related to cartilage. In vivo it is produced by chondrocytes and stored bound to the ECM. In unbound form it stimulates chondrocytes to synthesise ECM components such as collagen type II and proteoglycans (particularly aggrecan). It drives MSCs to differentiate into chondrocytes. The biochemistry of TGF-β has several subtleties: it is secreted in active or latent form as part of a 'large latent complex', whose active component is bound to a protein and a peptide and unable to bind to cell receptors until it is released from this complex. This arrangement enables the latent TGF-β to bind to and be stored in the ECM, for release only when chondrocyte stimulation is required (in response to ECM damage for instance). TGF-β activation typically occurs in response to mechanical or chemical stimulation [9,10]. Once activated it has a short half-life of only a few minutes as compared to the time-scale of days over which cell differentiation and proliferation occur [11]. Consequently it must bind rapidly to cells after activation to have any effect. For in vitro studies, TGF-β is replenished at regular intervals to ensure a sustained effect. For tissue engineering applications there are two dominant sources of TGF-β: it may be secreted by chondrocytes or it may be added exogenously. The key modelling question that we address in this paper is then to understand how these two sources interact and, in particular, to establish whether exogenous TGF-β alone is sufficient to drive stem cell differentiation.
A number of experimental strategies have been developed to drive the chondrogenesis of an in vitro population of MSCs. One strategy involves adding exogenous TGF-β to the initial population of stem cells [8]. Although TGF-β has a short half life, if its initial concentration is high enough it is possible that it will drive some cells to differentiate before it degrades completely. These newly differentiated chondrocytes will constitutively produce more TGF-β, and if they do so in a large enough quantity then this may drive further differentiation. Where too few cells are initially differentiated more TGF-β may be added to achieve complete differentiation. Mathematical modelling of this excitable system will help to identify conditions under which the system can be driven to constitutively produce enough TGF-β to sustain cell differentiation after the initial source of exogenous TGF-β has been depleted.
An alternative strategy for obtaining a population of differentiated cells is to co-culture MSCs with some harvested chondrocytes. There is evidence that co-cultures of this type produce cell populations with a more stable phenotype [12]. Under these conditions, the mechanism that induces chondrogenesis is similar to that above, except that the only source of the TGF-β is that produced by the chondrocytes. possible that even fewer chondrocytes and a lower concentration of exogenous TGF-β will be required than if the these strategies are used in isolation [13]. Mathematical models of this scenario can be used to assess the relative importance of these alternative mechanisms. Ideally experiments of the type shown in Figure   2 could be used to validate the model, but at present the available data consists of histology images at two sparse time points (at the start and end of the experiments), whereas model validation and parameter estimation would require quantitative data at multiple time points. This is beyond the scope of the present study but could form the basis for future work in which the mathematical model here is used to identify the time points at which the data is collected.
The purpose of this paper is to develop a mathematical model that describes the impact of TGF-β and chondrocytes on MSC proliferation and differentiation under various experimental conditions [8,12,13], and in doing so to identify the contributions from the various cellular and biochemical processes involved.
A qualitative comparison is appropriate since the reported results in experiments are largely qualitative in nature, typically based on stainings for ECM components of the type shown in Fig. 2, and therefore not suitable when making quantitative comparisons. In formulating our model we adopt a highly simplistic view of the relevant biology, assuming well-mixed populations of cells and spatially uniform chemical concentrations, as well as simplifying aspects of the TGF-β life cycle, while retaining all the key mechanisms involved in chondrogenesis. Previous models have considered TGF-β decay and activation in the context of fibroblast differentiation [14], with a detailed focus on receptor kinetics, as well as modelled in detail the TGF-β signalling pathway [15]. We do not attempt to include this level of detail in our model, instead focussing on the interactions between extracellular TGF-β and the two cell types.
The remainder of this paper is organised as follows. In Section 2, we develop our mathematical model and present some motivating simulations. The chondrogenesis strategy of co-culturing MSCs and chondrocytes is investigated in Section 3, where we introduce a fast kinetics approximation and derive an analytic expression for the critical density of chondrocytes required to drive chondrogenesis. The effect of adding exogenous TGF-β on a population of cells is explored by considering early-time behaviour in Section 4, and we obtain the critical concentration of added TGF-β required to drive chondrogenesis with a combination of analytic and numerical techniques. Finally, a hybrid strategy where exogenous TGF-β is added to a co-culture is investigated in Section 5 with a combination of analytic and numerical techniques. We discuss these results and present our conclusions in Section 6.

Model and biological background
We develop a time-dependent mathematical model for TGF-β mediated proliferation and differentiation of MSCs to chondrocytes. We neglect cell motility and any spatial variation in the densities. We assume that chondrocytes do not proliferate or dedifferentiate, which is a reasonable assumption for tissue engineering applications, where the chondrocytes are typically derived from articular cartilage and cultured in vitro [16,17,18]. In this context it is also reasonable to neglect cell death of both MSCs and chondrocytes on the relatively short time scale of interest (approximately 1-2 weeks) [18].
We propose an ordinary differential equation model to investigate how the cytokine, TGF-β, can drive chondrogenesis in a well-mixed (spatially uniform) population of cells; the dependent variables of this model are summarised in Table 1 chondrocytes: ologically appropriate description of the cell response to receptor bound TGF-β; following [19], for instance, the choice of Heaviside functions here has the advantage that it will permit further mathematical analysis.
The constants Λ 1 and Λ 2 are the rates at which differentiation and proliferation, respectively, occur where these effects are active. Here we assume that f d < f p , so that proliferation only occurs if the MSCs are also differentiating; this is broadly in line with the observed behaviour of MSCs that have been cultured in vitro, in a hydrogel, for a period of 1-2 weeks in the presence of the isoform TGF-β1 [5,18,20,21,22] in tissue engineering studies. We note that in other contexts and for other isoforms of TGF-β the effect may be different; for instance TGF-β in increases the proliferation of MSCs in bone formation [23], but has no effect on the proliferation of genetically engineered stem cells [24]. We further assume that Λ 1 > Λ 2 , so that when both effects are active, proliferation does not outpace differentiation. This is in line with observations of the relative rates of these effects in experimental studies [20,18]. We note that when Λ 1 < Λ 2 the model is physically unrealistic as the MSC population grows without bound.
In practice TGF-β can exist in up to eight different forms [9]. For simplicity, we distinguish just three forms of TGF-β and also track the quantity bound to the MSCs and chondrocytes. TGF-β is secreted as chondrocyte unbound latent TGF-β, c(t) unbound active TGF-β, a(t) (4b) (4a) a large, latent molecular complex, which comprises the TGF-β ligand, a peptide which prevents interaction with cell receptors, and a protein which allows the complex to bind to the ECM. For our purposes latent TGF-β can exist in two states; an unbound form, which we denote by c(t) (mass/volume), and an ECMbound form, which we denote by b(t) (mass/volume). In order for the ligand to be accessible to the cell receptors it must be activated, that is cleaved from the latent complex: we denote the concentration of active TGF-β by a(t) (mass/volume). A schematic diagram illustrating this life cycle is depicted in Figure 3.
Activation can arise from mechanical interaction with cell integrins [10] or through chemical interactions with proteases, thrombospondin-1, reactive oxygen species or extremes in pH [9]. Here, for simplicity, we assume that activation is caused by chemical interactions. Active TGF-β has a short half-life and so must bind rapidly to cell receptors or it will be degraded. As above, the fraction of bound TGF-β receptors per MSC for all cells in the population is denoted as f (t), and the fraction of bound receptors per chondrocyte for all cells in the population is denoted by g(t). Here, receptor-ligand dynamics are much simplified, as compared to other models [14,15], since our interest is in cell differentiation rather than the details of these dynamics. We assume that when an MSC differentiates the bound TGF-β remains bound and, therefore, MSC differentiation will give rise to a loss of MSC-bound TGF-beta and an equal and opposite source of chondrocyte-bound TGF-beta. We further assume that each MSC or chondrocyte has a fixed number of receptors; so that if all receptors are occupied the total bound mass of TGF-β per MSC is denoted F tot , and the equivalent quantity for the chondrocytes is denoted G tot . The receptor bound TGF-β for both cell types is internalised at a constant rate, and we assume that when a ligand-receptor complex is internalised it is replaced by a free receptor.
As for the cell densities, we assume that the system is well-mixed so that the time evolution of the various

Dependent variable Description
Units m density of mesenchymal stem cells (MSCs) Mio/mL n density of chondrocytes Mio/mL c concentration of unbound, latent TGF-β ng/mL b concentration of bound, latent TGF-β ng/mL a concentration of unbound, active TGF-β ng/mL f fraction of bound TGF-β receptors per MSC g fraction of bound TGF-β receptors per chondrocyte - Table 1: Description of the model variables, along with their units. Here "Mio" is an abbreviation of "million".
forms of TGF-β can be described by time-dependent ordinary differential equations (ODEs). The governing equations, along with a description of each term in underbraces, are presented in turn below: bound TGF-β to MSCs: loss of bound TGF-β due to MSC differentiation (6) bound TGF-β to chondrocytes: gain of bound TGF-β due to MSC differentiation (7) The dependent variables in Equations (1)-(7) are listed in Table 1 while descriptions of the model parameters, along with typical values, are given in Table 2. Most of these parameter estimates are taken directly from a combination of experimental studies and previous models of systems involving TGF-β [14,15]. Others are estimated from experimental studies in combination with values from other studies. For example, the estimate of λ 3 , is obtained from data given in [16] using parameter values from [25]. Similarly, in the absence of data with which to estimate the rate of differentiation, Λ 1 the value in Table 2 was chosen to be consistent with the time scales over which in vitro studies of these phenomena are typically conducted. Estimates are not given for the proliferation rate Λ 2 or the threshold values, f d and f p ; determining realistic values of these parameters is an aim of the present study. Note that the 'binding to cell receptor' terms in Equation  (5) are multiplied by F tot and G tot , whilst the corresponding terms in Equations (6) and (7) are not. This is because the units of f m, gn and a are different, and to ensure that all the terms in Equation (5) where stars denote dimensionless quantities. Here a 0 and m 0 are typical values for the concentration of TGF-β and cell densities; reasonable values for these scales, in the context of tissue engineering [3], are a 0 = 1 ng/mL and m 0 = 1 Mio/mL. The timescale of interest Λ −1 1 is that associated with cell differentiation; we shall see that the TGF-β dynamics operate on a much faster timescale than the timescale for differentiation.
Under these scalings, Equations (1)- (7) become where we have dropped stars, and expressions for the dimensionless parameters are given in Table 3, along with representative values. Note that many of these parameters are much greater than 1, indicating that the dynamics of the transfer between the various forms of TGF-β, or their decay/internalisation, are much faster than MSC proliferation and/or differentiation. We exploit this separation of time scales between the chemical and cell processes by introducing the inverse of the decay rate of the active form as small parameter, = 1/λ 8 , and scale the remaining large parameters as Substituting the rescaled parameters into (9)-(15) we obtain We note that the time derivatives in (19)-(23) are multiplied by the small parameter , reflecting that these dynamics occur on a fast time scale, while the cell processes represented by (17)-(18) act on a slower time scale.

Simulations of the full model
In the analysis that follows we shall see that the separation of time-scales between the chemical and cell processes leads to a rich solution structure. We begin by presenting numerical simulations of the full system The final yield of chondrocytes, n ∞ , depends on the initial density of seeded chondrocytes n(0) as shown in Fig. 4 and in more detail in Fig. 5(a). If the initial density is less than some critical value, n crit , then no MSC differentiation occurs, and so the final yield of chondrocytes is simply the initial density (n(0) = n ∞ ).
If the initial density exceeds n crit then eventually all the MSCs differentiate, so that n ∞ ≥ m(0) + n(0).
As shown in Fig. 4(a), as n(0) is increased beyond n ∞ there is a commensurate increase in the final yield n ∞ due both to the extra initial seeded chondrocytes and a small amount of MSC proliferation. Practical interest lies in determining the value of n crit , since there is no advantage in seeding with more chondrocytes than is necessary to get the system into a state where all MSC will differentiate, as well as in understanding how n crit relates to other parameters in the model. For example, as shown in Fig. 5(b) the value of n crit increases as the threshold value for differentiation f d is increased.

Strategy 2: Exogenous TGF-β driving MSC differentiation
We now consider the case where exogenous TGF-β is added to a population of MSCs. Here, the initial  Table 3. When the initial concentration of TGF-β is increased (from a(0) = 1 to a(0) = 5; dashed curves in Figure   6) more cells differentiate in the short time before the active TGF-β decays compared to the case where a(0) = 1. The larger chondrocyte population constitutively produces enough TGF-β (f > f d ) to ensure all MSCs eventually differentiate. For completeness, we note that if a very low initial concentration of a is used, or f d is large, then it is possible that no differentiation will occur.
The final chondrocyte yield n ∞ is shown for a range of initial concentrations a(0) in Fig. 7(a), showing three possible outcomes. If the initial concentration is less than some critical value a crit,1 , then the TGF-β degrades before any MSC differentiation takes place, giving a long time cell population consisting of the initial MSCs and no chondrocytes, and all forms of TGF-β zero. For moderate size concentrations, where a crit,1 < a(0) < a crit,2 , a few cells will be differentiated at early times but these are insufficient to drive  Table 3; thus the value of n crit is the same as in Fig. 5.  Fig. 7(b)).

Strategy 1: Chondrogenesis in co-cultures of MSCs and chondrocytes
The simulations presented in Figures 4 and 5 reveal that a population of MSCs can fully differentiate when they are co-cultured with a small number of chondrocytes if the initial proportion of chondrocytes exceeds some critical value, n crit say. That is, for a given initial chondrocyte density n(0) two long term outcomes are possible: • No differentiation: n(0) < n crit ⇒ n ∞ = n(0); • Full differentiation: n(0) > n crit ⇒ n ∞ ≥ m(0) + n(0).
In experimental studies the percentage of chondrocytes in co-cultures typically ranges from 25-66% [12], although ideally experimentalists aim to use fewer chondrocytes, since they are difficult to harvest [13]. Note that the chondrogenesis experiment in Fig. 2 used only 20% chondrocytes and the numerical simulations in Section 2.1.1 predicted that even fewer chondrocytes were required to drive differentiation, in turn this suggests that (assuming our parameter values are reasonable) fewer chondrocytes may need to be harvested to successfully employ this strategy.
We now exploit the different time scales of the cell process and TGF-β dynamics to make a quasi-steady (or 'fast kinetics') approximation to (19)-(23). This will allow us to estimate n crit in terms of the model parameters.

Fast-kinetics approximation of the TGF-β dynamics
Since the TGF-β dynamics occur on a faster timescale than the cell processes it is appropriate to make a quasi-steady approximation to Equations (19)-(23) by taking the limit as → 0. Equations (19)- (23) then reduce to the following algebraic expressions relating the various components (where a subscript 'FK' denotes the quasi-steady, or 'fast-kinetics', approximation of a particular variable), , Equations (24)-(28) implicitly define the time-evolution of the variables c, b, a, f and g, via their dependence on m and n whose evolution is defined by Equations (17) and (18). From (24) and (25) we note that c F K (t) and b F K (t) are proportional to the chondrocyte density, n(t). The remaining variables a F K (t), f F K (t) and g F K (t) are interdependent; while no simple expression in terms of the cell densities exists, some further analysis is possible. Substituting expressions (25), (27) and (28) into (26) we obtain the following cubic equation for a F K (t) = a F K (m(t), n(t)): We note that for n > 0 this expression has one real positive root, since its discriminant is always positive, and two stationary points at which a F K < 0.
Following some straightforward manipulation of the fast kinetics approximation (24)-(28) we can determine n crit , the critical value for the initial concentration of chondrocytes above which differentiation of the whole population occurs. This is the value of n which gives f F K = f d , the threshold value of f above which differentiation occurs. This value of f F K is associated with a threshold value of a F K = a d , and from rearrangement of (27) we find Substituting this value into equation (28) we define g d =k n9 a d /(k n9 a d + k n10 ). By rearrangement of equation (26) we obtain n crit in terms of m(0) and the system parameters, namely The value of n crit in (31) is in excellent agreement with the value obtained from the numerical solution of the full governing equations (dashed curve in Fig. 5(b)). We note that (31) depends on all the system parameters given in Table 3, except for the proliferation rate Λ 2 . This analytic expression reveals the dependence of the critical value on the model parameters. For example, we note that an increase in k 6 , the rate at which the latent form of TGF-β is activated, leads to a commensurate decrease in n crit , reflecting that the constitutively produced TGF-β is being transferred more efficiently.

Strategy 2: Chondrogenesis induced by exogenous TGF-β
We now determine the conditions under which a population of MSCs will eventually be fully differentiated by adding exogenous TGF-β. Although any added TGF-β quickly decays, as seen in the numerical simulations in Figure 6, it may still have a significant effect on the system in this short time period. The route to chondrogenesis is that in this short time a few MSCs are differentiated, and the resulting chondrocytes drive the remaining MSCs ultimately to differentiate in the manner described in Section 3. The simulations in Section 2.1.2 revealed three possible outcomes arising from the addition of an initial concentration a(0) of active TGF-β to a population of pure MSCs: • No differentiation: a(0) < a crit,1 ⇒ n ∞ = 0; • Partial differentiation: a crit,1 < a(0) < a crit,2 ⇒ 0 < n ∞ < n crit ; • Full differentiation: a(0) > a crit,2 ⇒ n ∞ ≥ m(0).
We examine the behaviour of (17)-(23) at early times using perturbation series techniques and, in so doing, determine how many chondrocytes are produced in response to the exogenous TGF-β, which in turn allows us to approximate a crit,1 and a crit,2 for different values of the initial concentration a(0).

Early times
The short half-life of active TGF-β means that any MSC differentiation associated with exogenously added TGF-β occurs a short time after it has been introduced to the culture medium. To this end we now analyse Equations (17)-(23) at early times in response to the addition of a known amount of active TGF-β at t = 0. Our goal is to estimate how many chondrocytes (if any) are produced during this short time.
To determine whether differentiation and/or proliferation will occur we need to solve for f . Thus for this analysis it is convenient to use Equations (17) and (18) to rewrite Equations (22) and (23) as differential equations for the per cell quantities f and g rather than the quantity per volume f m and gn. These equations When active exogenous TGF-β is added to population of MSCs the initial conditions are In general, these initial conditions will be incompatible with the fast-kinetics approximation (24)- (28), and so there is boundary layer at t = 0. We scale into this layer by setting t = τ and rescale a with a parameter δ that indicates the size of its initial magnitude, so that a = δA with A ∼ O(1) and a 0 ∼ O(δ); in Section 4.1.1 we take δ = to consider small initial concentrations of a, and then in Section 4.1.2 take δ = 1 to examine moderate initial concentrations of a. Substituting these scalings into (17)-(21), (32) and (33) gives that is, at leading order, the cell densities do not change from their initial values. Continuing to O( ) gives Our primary interest is in whether the exogenous TGF-β induces differentiation of MSCs at early times and, if so, to establish how many chondrocytes are produced; this is determined by integration of (46) as where τ * is an O(1) parameter for the width of the boundary layer (an appropriate choice for this is τ * = 2 log(1/ ), namely the time at which an exponentially decaying initial concentration of a will have decayed by two order of magnitudes in ). Similarly Equations (37) and (38) supply the following differential equations at leading order: These leading order solutions are consistent with there being no chondrocytes at leading order, and therefore no constitutive production of latent TGF-β. The corresponding equations at O( ) are The leading-order (and first-order) forms of Equations (39)-(41) differ slightly depending on the choice of δ, the initial magnitude of a, and in the following analysis we consider the cases where the system is excited by different initial concentrations at δ = and δ = 1, which each reveal different aspects of the behaviour of the system.

Small initial concentration
Here we assume the initial concentration is small, so that a(0) = A(0) = α 0 , where α 0 ∼ O(1). We seek to estimate a crit,1 , the critical initial value of concentration a(0) below which no MSC differentiation occurs. To this end we assume that for these small initial concentrations no differentiation or proliferation occurs, that is f < f d and f < f p , which implies that n 1 = 0 and m 1 = 0 from (45) and (46). By a similar argument to that applied to solve Equations (48) and (49), it follows from (50) and (51) that c 1 = b 1 = 0.
The leading order equations for f and A are then Continuing to O( ) in f we have These solutions are in excellent agreement with the numerical solution of the full system. See, for example, It is straightforward to show that f 1 attains its maximum at some time t * max . We seek the value of the initial concentration a(0) = α 0 at which this maximum corresponds to the critical value for differentiation (ie. f 1 = f d ). Rearranging (54), we deduce that no differentiation will occur if This inequality provides an asymptotic estimate for the first of the critical values in Figure 7(a), that is the value that separates regimes in which there is no cell differentiation from regimes in which there is partial differentiation. In Figure 8(c) we show that this approximation is in excellent agreement (less than 6.22% relative error where a(0) < 1) with the numerical results generated from the full model when f d is small, and (as might be expected) underestimates the critical value in the limit as f d → 1. Equation (55) reveals that this critical value depends only on f d and the system parameters involved in the binding and internalisation of TGF-β to MSCs. This reveals, for instance, that a crit,1 increases with the internalisation rate k m10 , indicating that more exogenous TGF-β must be added to achieve any cell differentiation if TGF-β is taken up by the cells at a higher rate.
To summarise, this analysis confirms the numerical result from Section 2. When a moderate concentration of exogenous TGF-β is added to a culture medium there is typically some differentiation of cells at early times, as indicated in the simulations in Section 2.1.2. The analysis of the previous section provided an excellent approximation of a crit,1 , and we now aim to estimate a crit,2 , the critical concentration above which sufficient cells are produced to drive differentiation of the whole population. The simulations presented in Section 2.1.2 suggest that a crit,2 is approximately an order of magnitude larger than a crit,1 . Therefore, in what follows we will assume that a crit,2 = O(1).
We assume that the a(0) = O(1) and δ = 1. Under these assumptions, it follows that at leading order Equations (39) and (40) supply   Figures 8(a) and (b), respectively. The horizontal dashed line in Figures 8(b) is the threshold value f d = 10 −2 . This solution for f reaches, but does not cross this threshold, indicating that the value of a(0) = a crit . The approximation to this critical value from (55) is shown as a dashed line in Figure 8(c) against the numerical critical values (solid lines) and is nearly indistinguishable for small a(0) and f d 1.
Here, fp = f d + 4 × 10 −2 and all other parameter values are given in Table 3; thus the value of n crit is the same as in Fig. 5.  Table 3.
since both b and n are zero to leading order, and this system has no closed form solutions. It is possible, however, to bound the solution to A 0 by using the fact that 0 ≤ f 0 ≤ F 0 in (56), where F 0 =k m9 α 0 /(k m9 α 0 + k m10 ) from the fast-kinetics expression (27); that is, as the initial concentration of a is decaying, we assume that f does not exceed the quasi-steady value it would take when the a concentration of a = α 0 is present in the system. The upper and lower bounds on A 0 are obtained by setting f 0 = F 0 and f 0 = 0, respectively, in (56) which gives and the lower bound in particular is an excellent approximation to the numerical solution of a; see, for instance, Fig. 9(a) where the relative error between the approximation (59) and the numerical solution is less than 2% for t < 0.05. These bounds are then used in (57) to determine upper and lower bounds on f 0 .
To find the upper bound on f 0 we set A 0 = A U 0 in the first term on the right hand side (so that the effect of this term is maximised) of (57) and A 0 = 0 in the second term (so that the effect of this term is minimised) To obtain a lower bound we set the A 0 = A L 0 in the first term and A 0 = α 0 in the second term, thereby minimising the growth and maximising the decay of f 0 , so that (57) becomes These bounds for f 0 are compared with the numerical solution in Fig. 9(b) with the upper bound again being a particularly close approximation to the numerics with a relative error of at most 8.61% for t < 0.1.
Having bounded f 0 we can estimate the number of chondrocytes produced in response to a given initial  Table 3. estimated by finding the roots of f U 0 − f d = 0, and noting that two such roots exist provided a(0) > a crit,1 from (55) (that is, when some MSCs have been differentiated). Having obtained these two values, τ U a and τ U b , it follows from (47) that the number of differentiated cells is A similar procedure can be performed with the lower bound on f 0 to obtain a lower estimate n L 1 for the number of chondrocytes produced. Both these estimates are are in good agreement with the numerical solution at moderate amplitudes of a(0) (see Fig. 9(c)), and only diverge from the computed value at larger concentrations.
Finally, upper and lower estimates of the critical value a crit,2 correspond to the values of a(0) such that n U 1 (a U crit,2 ) − n crit = 0 and n L 1 (a L crit,2 ) − n crit = 0, respectively. It is straightforward to find these values numerically in MATLAB, and the value of a U crit,2 from the upper approximation (vertical dotted line in Fig. 9(c)) slightly underestimates the numerically obtained value of this parameter. In the example shown no approximation of a L crit,2 is possible since n L 1 does not reach the critical value n crit . The upper and lower approximations of a crit,2 are shown against the differentiation threshold parameter f d in Fig. 10. Both are excellent approximations to the numerical solution for moderate amplitudes of a(0), and a U crit,2 remains close to the numerical value even for larger amplitudes, with a maximum relative error in the upper bound of 15.5% for the largest amplitude of a(0) shown here. As in Fig. 9(c), no lower approximation a L crit,2 exists beyond a moderate value of f d as the lower bound on n 1 never reaches n crit , even for large values of a(0).
To summarise, this analysis replicates the behaviour seen in the numerical simulations from Section 2.1.2 where full differentiation of a MSC population is achieved if sufficient exogenous TGF-β is added. The key insight that this analysis yields is an approximate expression for the concentration above which full differentiation takes place, as well as revealing the dependence of this value on model parameters. This estimate has been calculated by determining how many cells are produced by an initial concentration of exogenous TGF-β as described in Equation (64).

Hybrid strategy: add exogenous TGF-β to a co-culture
A final chondrogenesis strategy combines the strategies from section 3 and 4, as in the experimental study of [12] where a mixed cell population of 25% chondrocytes and 75% MSCs was cultured with various doses of TGF-β up to 10 ng/mL. The intent of this combined strategy is that is will be more cost-effective and robust than relying on driving chondrogenesis by either the use of co-cultures, or the addition of exogenous TGF-β alone; where growth factors are expensive and harvesting of chondrocytes is difficult it is likely that a hybrid approach will offer numerous advantages.
This situation is straightforward to implement in our model. We first obtain the fast kinetics values for all the TGF-β related quantities for a given initial co-culture population (24)- (28). Here, we assume that the number of seeded chondrocytes is less than n crit , so the co-culture does not completely differentiate without any extra stimuli. The effect of adding extra exogenous TGF-β to this system is then examined by a using fast kinetics values for c, b, f , g, m and n as initial conditions in the full numerical simulation, with the initial value for a, a(0) = a FK + a exogenous , perturbed from its fast kinetics value by the concentration of the exogenously added TGF-β.
The results of some example simulations are shown in Fig. 11, where this hybrid strategy is compared to the strategy of adding exogenous TGF-β to a population of pure MSCs. As shown in Fig. 11(a), the value above which any extra chondrocytes are produced a crit,1 is almost unchanged by seeding of the 3% chondrocytes, but the value for a crit,2 is decreased by around a factor of 2. This indicates that combining Here the value of "hybrid a crit,2 " for the case where some chondrocytes were in the initial cell population is significantly lower than for the case without any chondrocytes initially. (b) The dependence of a crit,2 for a range of initial chondrocyte densities n(0) for a fixed value of f d = 10 −2 . The value of the hybrid a crit,2 decreases markedly as the number of seeded chondrocytes is increased. Full differentiation occurs if a(0) and n(0) are chosen so that they lie above the solid curve.
Here, fp = 5 × 10 −2 and all other parameter values are given in Table 3.
these two strategies does indeed dramatically decrease the amount of TFG-β required to induce the full differentiation of the population.
Some further results are shown in Fig. 11(b) indicating how the cut-off value for full differentiation, a crit,2 , changes as the proportion of initially seeded chondrocytes and the concentration of exogenous TGF-β is increased. For a fixed value of the differentiation threshold parameter f d = 10 −2 the amount of exogenous TGF-β required to induce complete differentiation drops by more than an order of magnitude for a co-culture seeded with 10% chondrocytes, as compared to the case with purely MSCs. This prediction seems reasonable in light of the modelling of the individual strategies in isolation in Sections 3 and 4, and could result in cost-savings through the use of fewer chondrocytes and lower of TGF-β concentrations. Experimentally testing of this mixed strategy would be an interesting direction for future work.

Discussion and conclusions
We have developed a model to examine the differentiation of a population of MSCs into chondrocytes in response to TGF-β. There are two key points of novelty in this model. Firstly, the concentrations of the various forms of TGF-β are modelled explicitly, so that it is possible to track the storage of TGF-β in the ECM in a latent form for later activation. Secondly, the key aspects of the receptor kinetics involved are captured in a simple way, namely by the use of step functions to trigger differentiation or proliferation when the TGF-β bound to an MSC exceeds certain thresholds; this obviates the need to model the TGF-β signalling pathway explicitly. As discussed in Section 2, the use of step functions is standard practice in models of this type and has the advantage that the resulting model is amenable to detailed mathematical analysis. Possible issues with our model arise from the large number of parameters required to describe the complex biochemistry of the system. Although the values of these parameters were estimated from the literature, there is naturally a degree of uncertainty about the accuracy of these estimates. Similarly, the assumption that the cell populations and chemical concentrations are well-mixed means that any spatial effects (diffusion of chemicals for instance) are not considered here.
Two commonly used experimental chondrogenesis strategies were simulated. Our analysis of chondrogenesis driven by co-culturing MSCs and chondrocytes in Section 3 revealed that there is a critical initial density of chondrocytes n crit , given in (31), required to induce a population of MSCs to completely differentiate. The implication of this result is that if the initial cell density is above this critical value then the long term cell population will consist only of chondrocytes, but if the initial density is below the critical value then no extra chondrocytes will be produced and the long term cell population will be identical to the initial population. Thus, the optimal strategy suggested by the model is to add an initial density of chondrocytes which just exceeds the critical value; there is no real benefit to adding more chondrocytes that is sufficient to trigger chondrogenesis in terms of the ultimate yield or time taken to fully differentiate the MSCs, and so harvesting more chondrocytes than required to achieve this is unnecessary. The existence of such a critical value is broadly in line with behaviour seen in co-culture experiments. Under the conditions used in the in vitro study [29] the value of n crit is between 5% and 25%, which is consistent with the values in Fig. 5(b).
This suggests that the parameter values in Table 2 are not unrealistic and that a threshold parameter for differentiation f d in the range 10 −3 -10 −2 corresponds well with the results of in vitro studies.
Similarly, our analysis of the strategy of adding exogenous TGF-β to a population of MSCs revealed that in this excitable system the final yield of chondrocytes is dependent on the initial concentration of TGF-β, and that a population of MSCs will be completely differentiated if sufficient exogenous TGF-β is added.
No cells are produced if the initial concentration is below a cut-off value a crit,1 , and an analytic estimate of this value (see Equation (55)) was obtained. The entire MSC population can be driven to differentiation if the initial concentration exceeds a threshold value a crit,2 , and a procedure to estimate this cut-off value, revealing its dependence on various modelling parameters, is given in Section 4.1.2. Our estimate of this second critical concentration, which is of particular interest in the context of tissue engineering, is slightly lower than the 10 ng/mL TGF-β conventionally added to in vitro experiments, suggesting that a similar outcome could be achieved experimentally with a lower concentration than is currently used. It would be interesting to test this experimentally, but this is beyond the scope of the present study.
Future extensions to this model could include an investigation of the effect of mechanical loading, for instance by using a piston to apply a load to cells within a hydrogel construct, as this kind of stimulation has been demonstrated to promote chondrogenesis [8] and has been proposed as a route to TGF-β activation [10]. Mechanical activation of this kind could be added to model of this paper by, for instance, replacing the chemical activation term in Equations (4) and (5) with a mechanotransductive term that depends on the local stress, meaning that TGF-β stored in the ECM could be activated and made available to cells in response to an applied mechanical load. This would necessarily require the extension of the well-mixed model of this paper, to one that includes spatial effects. A model which includes spatially varying cell densities and chemical concentrations would also allow other interesting biological features of the system to be investigated. For instance, this might include capturing the competition betweeen the diffusion of TGF-β and its rapid decay, which may mean that, for instance, TGF-β added to one side of a culture medium may not evenly spread throughout the medium, possibly resulting in isolated regions of the medium where no TGF-β is present. It would also allow consideratioin of cell movement by random motility and chemotaxis, both of which are likely be important in the context of cartilage tissue engineering to achieve a zonated cell distribution characteristic of cartilage, similar to that shown in Figure 1.
Since the synthesis of ECM components during chondrogenesis will dynamically change the bulk mechanical properties of extracellular environment [2,30], further work is needed to couple the mechanical loading model with a model for the production of ECM components, in particular collagen type II and glycosaminoglycan. These extensions demand consideration of spatial effects, rather than the well mixed approach used in this paper, and developing such an approach will allow examination of further experimental questions of interest. For instance, in the context of tissue engineering it is necessary to understand the balance between the diffusion of exogenously added TGF-β through a medium and the rate at which it binds to that medium. Spatial dependence would also enable determination of the kind of stimuli required to drive some initial population of cells into a zonated final distribution of chondrocytes, as in natural articular cartilage.
This model elucidates a key component of the mechanism behind chondrogenesis and has revealed the excitable nature of this system. This is of particular interest in the context of tissue engineering and has great potential to refine experimental practice by predicting the concentration of TGF-β, or the number of chondrocytes, required to achieve full differentiation of a population of MSCs. This could, for instance, to cost savings via a reduction in the amount of TGF-β being added to a culture medium, or by harvesting fewer chondrocytes.