A stochastic model of corneal epithelium maintenance and recovery following perturbation

Abstract

Various biological studies suggest that the corneal epithelium is maintained by active stem cells located in the limbus, the so-called limbal epithelial stem cell hypothesis. While numerous mathematical models have been developed to describe corneal epithelium wound healing, only a few have explored the process of corneal epithelium homeostasis. In this paper we present a purposefully simple stochastic mathematical model based on a chemical master equation approach, with the aim of clarifying the main factors involved in the maintenance process. Model analysis provides a set of constraints on the numbers of stem cells, division rates, and the number of division cycles required to maintain a healthy corneal epithelium. In addition, our stochastic analysis reveals noise reduction as the epithelium approaches its homeostatic state, indicating robustness to noise. Finally, recovery is analysed in the context of perturbation scenarios.

Appendices

Appendix A Parameter estimations

1.1 A.1 Mouse

We first note that the corneal circumference of a mouse is $\sim 10{,}000\,\upmu \hbox {m}$ (Di Girolamo et al. 2015; Dorá et al. 2015a; Dorà et al. 2015b) and a typical basal cell diameter is $\sim 10\,\upmu \hbox {m}$ (Romano et al. 2003). If stem cells simply formed a one-cell thick ring, a total stem cell population of $\sim 1000$ cells could be accommodated along the corneal-limbal border. Note, however, that an estimated 250–300 are active (Dorá et al. 2015a; Dorà et al. 2015b) at homeostasis. To accommodate scenarios that can range from healthy to pathological, or eye sizes from larger to smaller, we assume the number of SCa in the limbus ranges between 100 and 1000.

Although the cornea is dome-shaped, for the purposes of the model we have assumed it is a hemisphere with a circumference of approximately 10, 000 $\upmu \hbox {m}$. Then, the radius of the corneal is $r_{corneal}=1592$$\upmu \hbox {m}$, from which the corneal area is $A_{corneal}=2 \pi r^{2}_{corneal}$$\upmu \hbox {m}^{2}$. Similarly, the average area occupied by a basal corneal cell (assuming that the cell is a disc in the 2D plane) is $A_{cell}=\pi r^{2}_{cell}$$\upmu \hbox {m}^{2}$, where $r_{cell}=5$$\upmu \hbox {m}$. Thus, an estimate of cells that can fit in the corneal epithelium is given by:

$$\begin{aligned} \frac{A_{corneal}}{A_{cell}}= \frac{2 \pi r^{2}_{corneal}}{\pi r^{2}_{cell}}=202{,}757 \end{aligned}$$

(30)

and to take into account not just the normal conditions, we can introduce the magnitude of $10^{5}$ as a guideline baseline value for the number of cells required to populate a small cornea.

For mouse we have a number of sources that provide indications of stem cell and TAC division rates. If it is assumed that mouse limbal epithelial SCs are equivalent to BrdU “label-retaining cells”, which include slow-cycling stem cells, it can be estimated that certain limbal epithelial SCs do not divide more often than once per two weeks ($\sim 14$ days). This calculation follows from detectable BrdU retention for at least 10 weeks (Douvaras et al. 2013), and that BrdU is probably diluted to undetectable levels after 4–5 cell divisions (Wilson et al. 2008). However, this is quite likely to provide an approximate lower bound for division rates, as it remains quite possible that certain SCs divide significantly more quickly and may not be detected by the label-retaining cell approach. As such, the mean SC cell cycle time may be considerably less than 2 weeks. Of course division rates are ultimately bounded by the minimum length of time needed to complete the cell cycle, which would be of the order of several hours to a day. Consequently, we take a range 6 h to 16 days for (active) stem cell doubling times.

Experimental studies on the TAC cell cycle in the peripheral corneal epithelium indicate that almost $50\%$ of basal corneal epithelial cells are in S-phase of the cell cycle, during a 24-h labelling period (Urbanowicz et al. 2011). This suggests a minimum cell doubling time of just over 2 days but it would be longer if certain TACs cycle more slowly. Similarly, an average mitotic rate of $37\%$ of basal layer cells per day can be derived for rats from the results reported by Bertalanffy and Lau (1962) and this suggests a minimum cell doubling time of about 2.7 days. (The original results showed that $14.5\%$ of all corneal epithelial cells divided per day and results for the mouse imply that about $38.8\%$ of mouse corneal epithelial cells are in the basal layer (Douvaras et al. 2013)). Other experiments on the TAC cell cycle in the peripheral corneal epithelium have estimated it as approximately as 72 h for the mouse (Lehrer et al. 1998). Overall the results show that the average doubling time for TACs is about once every 2–3 days but may be longer in the central corneal epithelium (Lehrer et al. 1998). While we centre on an average rate of 2 days, for our studies we again use a range of 6 h to 16 days to include scenarios under normal and abnormal conditions.

1.2 A.2 Human

Experimental data suggests that the average corneal diameter in human eye is $11.71\pm 0.42 \,\hbox {mm}$, (Rüfer et al. 2005) implying a corneal circumference $\sim 36.770\,\hbox {mm}^{2}$. In the absence of specific data, we consider an analogous case to the mouse and suppose the circumference corresponds to the corneal-limbal border. Assuming limbal corneal cells are $10 \,\upmu \hbox {m}$ in diameter, we estimate that there is a room for $\sim$ 3000–40000 limbal cells forming a one-cell thick ring; although (in contrast to the mouse case) some biological studies suggest that they are asymmetrically distributed (Wiley et al. 1991; Pellegrini et al. 1999; Shanmuganathan et al. 2007). If a similar fraction (to that of mouse) of this population is taken to be active, we estimate $\sim 1000$ active stem cells ($SC_{a}$) in the human limbus. Again, we consider an order of magnitude range about this value ($\sim$ 400–4000).

Using the same calculations adapted from the mouse case gives an order of magnitude of $10^{6}$ basal epithelial cells fitting in the human cornea.

1.3 A.3 Rat and Rabbit

To demonstrate variability across other species, we note that rat and rabbit corneas have average diameters of $5.5 \,\upmu \hbox {m}$ (Cabrera et al. 1999) and $14.375 \,\upmu \hbox {m}$ (Tsonis 2011) respectively. Straightforward calculations show that the circumferences will be $17{,}270 \,\upmu \hbox {m}$ and $45{,}138\,\upmu \hbox {m}$ respectively. Making the same assumptions as earlier, this would allow for a total of 1727 and 4513 stem cells and, if again approximately 1 / 4 are active, $\sim 450$ and $\sim 1200$ active stem cells for rat and rabbit respectively. Calculating an estimate for the total number of cells that can fit into the basal epithelium yields a magnitude $\sim 10^{5}$ for rat and $\sim 10^{6}$ for rabbit, the former the same magnitude as the mouse and the latter similar to the human eye.

For a rabbit corneal epithelium, experimental data on TAC doubling time suggests once every 18 h (3 / 4) (Castro-Muñozledo 1994). We are lacking such data for the rat eye. Nevertheless, the parameter spaces provided throughout the paper can give a rough estimate of the TAC generations required for the epithelium maintenance for each of rat and rabbit eye.

Appendix B Derivation of matrices included in Lyapunov equation

1.1 B.1 Jacobian matrix

The Jacobian matrix J can be derived from the stochastic mean system (5)–(7) obtained in Sect. 2.3. Matrix J of our n-ODEs system for the stochastic means of TACs is:

$$\begin{aligned} \begin{aligned} \mathbf {J} = \begin{bmatrix} -\beta&\quad 0&\quad 0&\quad 0&\quad \ldots&\quad 0 \\ 2\beta p_{T,T}+\beta p_{T,TD}&\quad -\beta&\quad 0&\quad \ldots&\quad \ldots&\quad 0 \\ 0&\quad 2\beta p_{T,T}+\beta p_{T,TD}&\quad -\beta&\quad 0&\quad \ldots&\quad 0\\ \vdots&\quad \ddots&\quad&\quad \ddots&\quad&\vdots \\ \vdots&\quad \vdots&\quad \ddots&\quad&\quad \ddots&\quad \vdots \\ 0&0&\ldots&\quad 2\beta p_{T,T}+\beta p_{T,TD}&\quad&-\beta \\ \end{bmatrix} \end{aligned} \end{aligned}$$

Note that $J_{ij}= \dfrac{\partial }{\partial \phi _{j}}(\partial _{t} \phi _{i})$ where $\phi _{i}=N_{T_{i}}$ and $\phi _{j}=N_{T_{j}}$ with $j=1,\ldots ,n$.

1.2 B.2 Stoichiometric matrix

For the stoichiometric matrix we are only interested in the number of TACs at each reaction. Denoting the reactions as $r_{k,l}$ with k the reacting population (i.e. $k=0,1,\ldots ,n$, where $k=0$ corresponds to the $SC_{a}$ division to $TAC_{1}$ and $k=i$ the $TAC_{i}$ divisions) and the pathway indicator is l (hence $l=1,2,3$). The reactions can be written as

$$\begin{aligned} r_{0}:&\emptyset \xrightarrow {\alpha } TAC_{1} \end{aligned}$$

(31)

$$\begin{aligned} r_{i,1}:&TAC_{i}\xrightarrow {\beta p_{TD,TD}}\emptyset \end{aligned}$$

(32)

$$\begin{aligned} r_{i,2}:&TAC_{i}\xrightarrow {\beta p_{T,T}}2TAC_{i+1} \end{aligned}$$

(33)

$$\begin{aligned} r_{i,3}:&TAC_{i}\xrightarrow {\beta p_{T,TD}}\emptyset + TAC_{i+1} \end{aligned}$$

(34)

$$\begin{aligned} r_{n}:&TAC_{n}\xrightarrow {\beta }\emptyset \end{aligned}$$

(35)

where $i=1,\ldots ,n$ denotes the number of TAC generation. The stoichiometric vector for $TAC_{1}$ is [1 −1 −1 −1 0 $\cdots$ 0], for $TAC_{i}$ is [0 $\cdots$ 0 2 1 −1 −1 −1 0 $\cdots$ 0] and for $TAC_{n}$ is [0 $\cdots$ 0 2 1 −1]. As an example for the stoichiometric matrix $\mathbf {S}$, let us assume that the total number of TAC generations is 3, then

$$\begin{aligned} \mathbf {S}= \begin{bmatrix} 1&\quad -1&\quad -1&\quad -1&\quad 0&\quad 0&\quad 0&\quad 0 \\ 0&\quad 0&\quad 2&\quad 1&\quad -1&\quad -1&\quad -1&\quad 0 \\ 0&\quad 0&\quad 0&\quad 0&\quad 0&\quad 2&\quad 1&\quad -1 \end{bmatrix} \end{aligned}$$

Note that in the stoichiometric matrix for n TAC generations, the number of zero elements at the start of each row (excluding the first and last row which correspond to the first and last TAC generation respectively) will be $N_{0}=3(i-1)+2$ with $i=2,\ldots ,n-1$. Hence, the position of the first non-zero element in each row ($i=2,\ldots ,n-1$ ) follows the sequence $\sum _{i=2}^{n-1}3(i-1)$.

1.3 B.3 Vector of macroscopic rates

To find the vector of macroscopic rates we recall the reactions 31-35 listed in Appendix 1 with corresponding rates:

$$\begin{aligned} (\alpha [SC_{a}], \beta p_{TD,TD}(i)[TAC_{i}], \beta p_{T,T}(i)[TAC_{i}], \beta p_{T,TD}(i)[TAC_{i}], \beta [TAC_{n}]).\nonumber \\ \end{aligned}$$

(36)

Hence, the vector for our system is

(37)

1.4 B.4 Diffusion matrix

For the elements of the diffusion matrix, as already discussed in the text, we used

$$\begin{aligned} \mathbf {D}=\mathbf {S} \cdot \mathbf {F} \cdot \mathbf {S^{T}}, \end{aligned}$$

(38)

inside the matlab code, with $\mathbf {F}$, $\mathbf {S}$ and $\mathbf {S^{T}}$ determined as above.