Spatiotemporal Dynamics of Insulitis in Human Type 1 Diabetes

Type 1 diabetes (T1D) is an auto-immune disease characterized by the selective destruction of the insulin secreting beta cells in the pancreas during an inflammatory phase known as insulitis. Patients with T1D are typically dependent on the administration of externally provided insulin in order to manage blood glucose levels. Whilst technological developments have significantly improved both the life expectancy and quality of life of these patients, an understanding of the mechanisms of the disease remains elusive. Animal models, such as the NOD mouse model, have been widely used to probe the process of insulitis, but there exist very few data from humans studied at disease onset. In this manuscript, we employ data from human pancreases collected close to the onset of T1D and propose a spatio-temporal computational model for the progression of insulitis in human T1D, with particular focus on the mechanisms underlying the development of insulitis in pancreatic islets. This framework allows us to investigate how the time-course of insulitis progression is affected by altering key parameters, such as the number of the CD20+ B cells present in the inflammatory infiltrate, which has recently been proposed to influence the aggressiveness of the disease. Through the analysis of repeated simulations of our stochastic model, which track the number of beta cells within an islet, we find that increased numbers of B cells in the peri-islet space lead to faster destruction of the beta cells. We also find that the balance between the degradation and repair of the basement membrane surrounding the islet is a critical component in governing the overall destruction rate of the beta cells and their remaining number. Our model provides a framework for continued and improved spatio-temporal modeling of human T1D.

1. Introduction

Type 1 diabetes (T1D) is an auto-immune disease characterized by the selective destruction of pancreatic beta cells in the islets of Langerhans by the immune system (Eisenbarth, 1986; Atkinson, 2012; Boitard, 2012; La Torre and Lernmark, 2012; Pugliese, 2014; Richardson et al., 2014; Roep and Tree, 2014). This destruction takes place during an inflammatory phase, known as insulitis, in which various immune cells infiltrate the islets (Lecompte, 1958; Gepts, 1965; Willcox et al., 2009; Morgan et al., 2014). As the beta cell mass decreases over the course of the disease, the ability of the islets to secrete sufficient quantities of insulin to properly regulate blood glucose levels becomes compromised. As a result, patients with type 1 diabetes ultimately become reliant on the lifelong administration of external insulin.

Well-defined genetic components have been identified which predispose individuals to T1D; in particular the HLA-genotype (Itoh et al., 1993; Somoza et al., 1994). However, this alone is not sufficient to predict which individuals will develop the disease with precision. Since monozygotic twin studies have shown limited pairwise concordance for T1D (Barnett et al., 1981; Lo et al., 1991; Redondo et al., 2001), it is clear that environmental factors, such as viral infection, vitamin D status and childhood nutrition may also contribute significantly to the development of T1D (Knip et al., 2005).

Among the different subtypes of immune cells which infiltrate islets, CD8+ T cells are considered as the likely mediators of beta cell destruction (Bottazzo et al., 1985; Itoh et al., 1993; Somoza et al., 1994). It is widely believed that these promote beta cell apoptosis by both direct and indirect mechanisms and that macrophages then clear dying and dead beta cells very quickly. Other T lymphocytes, such as those expressing CD4, are also thought to play a role, though their precise functions are less clear (Willcox et al., 2009; Richardson et al., 2011). In addition to the T cells, B lymphocytes (CD20+) are also present in significant numbers during certain stages of insulitis in some patients. Indeed, recent evidence has suggested that the number of B cells and/or the ratio of B-cells to CD4+ cells present in the infiltrate can be used to classify the disease into two distinct phenotypes. These have been defined as “hyper-immune,” characterized by elevated numbers of CD20+ cells and a rapid loss of beta cell mass whereas, by contrast, the “pauci-immune” phenotype, is associated with a lower proportion of B cells and a much slower destruction of beta cells (Morgan et al., 2014; Leete et al., 2016). The mechanisms by which the B cells affect the rate of disease progression are unknown, but it is possible that they collaborate with CD8+ T cells to drive beta cell loss (Huppa and Davis, 2003).

Progress in understanding the cellular and molecular mechanisms underlying insulitis in humans has been hindered by the paucity of available samples from patients who died at, or close to, disease onset. Fewer than 200 such samples are available worldwide (Gepts, 1965; Foulis and Stewart, 1984; Klöppel et al., 1985; Dotta et al., 2007; Walker et al., 2011; Campbell-Thompson et al., 2012; Pugliese, 2014), and inferring the time course of a disease process from histological samples is fraught with difficulty since a range of assumptions and extrapolations about the likely progression are inevitably required to achieve this.

To offset this problem, much work has been performed using animal models such as the non-obese diabetic (NOD) mouse as a proxy for the human condition (Kachapati et al., 2012). The advantages of this are clear—experiments can be performed over relatively short periods of time and analysis of circulating and pancreatic lymphyocyte populations is achieved more readily during the course of disease. Such studies have, therefore, been extremely informative as a means to identify important cellular and molecular factors involved (Lally and Bone, 2003), although considerable effort needs to be made to verify that results are translatable to humans. While a large number of potential therapies have been identified in rodents, a means to prevent the human disease remains elusive (Brehm et al., 2012; In't Veld, 2014; Pugliese et al., 2014; Reed and Herold, 2015).

An alternative approach, which may shed light on human insulitis is offered by mathematical modeling. By constructing environments that mimic the pancreas and immune system, experiments can be performed in silico, and the results compared to empirical data arising from studies of relevant human tissue. Through model construction and analysis, candidate mechanisms giving rise to T1D can be interrogated in a systematic way. To this end, a number of authors have developed models describing: beta cell function (Bertram and Sherman, 2004), defective macrophage clearing (Marée et al., 2006), immune cell populations (Mahaffy and Edelstein-Keshet, 2007), multi-clonal populations of immune cells (Khadra et al., 2009, 2011; Jaberi-Douraki et al., 2015), immune cell cycles (Jaberi-Douraki et al., 2014b), and apoptotic stress generated by the loss of beta cell mass (Jaberi-Douraki et al., 2014c).

A common theme among the models outlined above is that they can be regarded as “lumped models,” in that they deal with averaged quantities taken over the whole body. For certain applications, this seems appropriate, as it is commensurate with current clinical practice. The biomarkers clinicians have available are typically derived from blood samples, in particular, measures of HbA1c and fasting glucose levels (American Diabetes Association, 2006; Inzucchi, 2012). Other biomarkers, such as measures of C-peptide (NIC, 2015) and islet cell and anti-insulin antibodies (Taplin and Barker, 2008) can be used to identify individuals who may be susceptible to diabetes and, though it is not currently routine clinical practice, they can also be used to aid classification of diabetes (Jones and Hattersley, 2013). Biomarkers obtained from blood samples are, by their nature, whole body measures. However, it is known from histological samples that insulitis displays a pronounced spatial dependence (Willcox et al., 2009; Morgan et al., 2014). This is due in part to the distribution of the islets within the pancreas, but is also a consequence of the presence of a basement membrane around the islets (Korpos et al., 2013), establishing a barrier to immune cell infiltration. Interactions and communication between immune cells also play a role and are likely to contribute to this spatial dependence.

Insulitis is a spatially heterogeneous process, both within an individual islet and across the pancreas as a whole. Islets that are heavily infiltrated can be located near to islets that are free from inflammation. Moreover, within an infiltrated islet, the destruction of beta cells appears not to follow a “wave-like” profile, as might be expected, but seems more random (Willcox et al., 2009; Morgan et al., 2014). These observations lend credence to the notion that immune cell communication is a critical component of insulitis.

To investigate insulitis and its relevance to the progression of T1D, these spatial aspects must thus be taken into account. Moreover, since the number of immune cells infiltrating any given islet is low (Willcox et al., 2009; Morgan et al., 2014), density based approaches, such as ordinary differential equation (ODE) modeling, cannot be applied with precision and alternatives are required.

In this manuscript, we construct an agent-based model of the insulitis process within a single islet. The model is posed on a regular domain that represents the space surrounding the islet. The principal agents in our model are the immune cells and we describe how they locate and target beta cells and how they interact with one another. By using an agent-based approach, we avoid the problems associated with density-based (averaged) descriptions of cell mass. We begin by introducing the relevant biology for our problem. Following this, we describe the development of, and results from, the agent-based model simulations. Finally, we conclude with a discussion of the model and its potential extensions.

2. Relevant Biology

In this section, we highlight the biology corresponding to the components of our model.

2.1. Cell Types

Our model consists of three distinct cell types—two types of lymphocyte and the beta cells. We consider in particular CD8+ T lymphocytes and CD20+ B lymphocytes. Of the various immune cells that have been shown to be involved in human insulitis, these are the ones that are believed to be most relevant to the questions under investigation in our study. The CD8+ cells are implicated in the destruction of the beta cells, while the CD20+ cells are deduced to play an important role, since their absence is associated with a weakly aggressive phenotype (Morgan et al., 2014). In our model, the B cells will essentially act as antigen presenting cells to the T cells (Rodríguez-Pinto, 2005). We assume that macrophages efficiently clear apoptotic beta cells. We do not account for mitosis in any of the cell types; whilst division and selection of relevant lymphocytes has been demonstrated in the lymphatic system, evidence for it in the peri-islet space is minimal (Willcox et al., 2010). Immune cells are assumed to have a finite lifespan and beta cells die following interactions with T cells (Cnop et al., 2010). We do not consider other cell types, such as α and δ cells, that are also found within the islets of Langerhans (Kim et al., 2009).

2.2. Basement Membrane

Individual islets are contained within an encapsulating membrane. This membrane is comprised of various laminins and collagen (Korpos et al., 2013) and acts as barrier to invading immune cells. In mice with T1D, these compounds are lost, suggesting that the basement membrane has been degraded. In our model, we assume that the T cells are responsible for degrading the membrane, but in principle, other cells could be also responsible for this. We also make the assumption that the membrane can be repaired over time.

2.3. Cell Movement

We assume that beta cells have a fixed location, but that immune cells are free to move around. The direction of this movement is mediated by a chemokine gradient, which will be described below (Stein and Nombela-Arrieta, 2005). In addition, immune cells of different types are attracted to one another, instigated by different chemical signaling pathways (Janeway et al., 1985). Cells are not allowed to pass through the basement membrane, nor through each other, but are allowed to remain in contact with one another.

2.4. Chemokine

The beta cells are considered to secrete a chemokine signal that attracts immune cells toward the islet (Christen and Von Herrath, 2004). In our study, we assume that the chemokine molecules are bound to the membrane of the beta cells and may also be cleaved off so that they diffuse freely.

2.5. Cell Interactions

In their role as antigen presenting cells, the B cells can cause the T cells to enter an activated state, in which both their sensitivity to the chemokine gradient and their killing efficiency are enhanced (Friedl et al., 2005). We assume that activated T cells have a shorter lifespan than their non-activated counterparts (Green et al., 2003).

When a T cell encounters a beta cell, it triggers a pro-apoptotic pathway within the beta cell, ultimately resulting in the death of the beta cell (Cnop et al., 2005).

The behavior of the cells and membrane are summarized below.

2.6. Summary of Basic Behavior

Below, we will summarize the behavior of each of the cells types included in our modeling.

T cells

• Move up chemokine gradients

• Degrade basement membrane

• Kill beta cells

• Interact with B cells

When activated by B cells, T cells become more sensitive to the chemokine gradient and become more effective at killing beta cells.

B cells

• Move up chemokine gradients

• Form productive interactions with T cells

• Activate T cells

Beta cells

• Immobile

• Killed by T cells

We do not include in our model the possibility of beta cell replication. Studies have suggested that beta cell proliferation is increased, particularly in the early stages of insulitis (Willcox et al., 2010). Here, we disregard this effect.

Basement membrane

• Degraded by T cells.

• Can self-repair.

3. Methods

3.1. Building a Cellular Automaton

ODE-based modeling approaches average the dynamics associated with insulitis over space essentially representing the entire pancreas. Such approaches have provided numerous insights into the potential roles of T and B cell involvement in T1D and we refer the reader to Jaberi-Douraki et al. (2014a) for a recent and comprehensive review of this literature. However, we know that insulitis is heterogeneous and that the number of immune cells directly involved in the destruction of the beta cells in humans is low - far lower than would be necessary for a density approximation to be justified (Willcox et al., 2009; Morgan et al., 2014). Instead, we employ a different, agent-based, modeling approach (An et al., 2009; Holcombe et al., 2012). Our agents, the immune and beta cells, occupy finite area within a two dimensional space. In a real islet, cells are free to move in three dimensional space, but for ease of modeling (as a starting point of first approximation) and for a more direct comparison with data, we shall here restrict ourselves to the planar case. The agents are given specific rules to locomote, interact with one another and with the basement membrane. Simulations will then be performed to track emergent behavior arising from the basic rules. In this way, the model will evolve over a spatially extended domain in which cell-cell interactions can be explicitly defined.

Note that the rules that govern our agents' dynamics are not intended to be exhaustive lists for all potential behaviors exhibited by the cells. Instead, they are intended to represent a minimal set of interactions that can explain the observed phenomena. Additionally, the model can be iteratively developed to account for new interactions as required.

3.2. Cell Movement

In this modeling framework, the cells are treated as being discrete circular particles (Levine et al., 2000; Palsson and Othmer, 2000; Tijskens et al., 2003; Maini and Baker, 2011; Bruna and Chapman, 2012; Plank and Simpson, 2012). Cell movement is simulated by modeling the forces that act upon cells from different sources and then by resolving these forces using Newton's second law of motion. In the absence of external forces, intrinsic forces acting upon the immune cells cause them to diffuse randomly, obeying a random walk. In the presence of a chemokine gradient, these same cells will seek to move up the gradient. Cells may be attracted to one another and may be in contact with one another for extended periods of time, that is, we allow cells to overlap, but they may not pass through each other.

3.3. Immune Cell Forces

Before discussing the governing equations, it is useful to define indexing sets for each of the cell types. In the following, we define the sets ${T},{B}$, and β as sets containing the indices for T, B, and beta cells respectively. From Newton's second law of motion, we have:

$m_i\frac{{\text{d}}^2{x}_i}{\text{d}{t}^2}=-{\eta }_i\frac{\text{d}{x}_i}{\text{d}t}+F_i+\sigma {\xi }_i(t),i\in \mathcal{T}\cup ℬ,(1)$

for each cell where $x_i\in {ℝ}^2$ is the cell's location, m_i ∈ ℝ_≥0 is the mass of the cell, η_i ∈ ℝ_≥0 is the cell's viscosity and $F_i\in {ℝ}^2$ is the force acting on the cell. The final term in Equation (1) represents a Gaussian white noise process: 〈ξ(t)〉 = 0, 〈ξ(t)ξ(t + s)〉 = δ(t – s), with strength σ such that the immune cells perform random walks in the absence of any other forces (Palsson and Othmer, 2000; Middleton et al., 2014). In biological systems, we assume that the cells have low Reynold's numbers such that inertial forces are small compared to viscous ones. Under this approximation, we can replace (Equation 1) with:

$\frac{\text{d}{x}_i}{\text{d}t}=\frac{1}{{\eta }_i}{F}_i+\sigma {\xi }_i(t),i\in \mathcal{T}\cup ℬ$

For simplicity, we set η_i = 1 for all cells. Note that, in a general mathematical framework, specific values for η_i can be absorbed into the definition of F_i. For our model definition, we include heterogeneity between immune cells through variations in their sensitivity to chemical gradients, as we shall discuss in Section 3.11. We will also initially assume that the cells within the subgroups are homogeneous with respect to their geometry, that is, they share a common radius r_imm.

The force F can be broken up into its constituent parts:

$F_i=F_i^{\text{chemo}}+F_i^{\text{cell-cell}}+F_i^{\text{islet}}+F_i^{\beta }.(2)$

In the above, $F_i^{\mathrm{\text{chemo}}}$ represents chemotactic force, $F_i^{\mathrm{\text{cell-cell}}}$ represents cell attraction and repulsion, whilst $F_i^{\mathrm{\text{islet}}}$ represents the interaction with the basement membrane and $F_i^{\beta }$ represents interaction with the beta cells.

3.4. Chemokine Signaling

Since the chemokine is a chemical signal, we establish a gradient using a reaction-diffusion equation:

$\frac{\partial C}{\partial t}=D{\nabla }^{2}C+f(C).(3)$

We assume that the chemokine evolves on a much faster timescale compared to the movement of the cells, so we assume the chemokine to be at steady state by setting the LHS of Equation (3) to zero. The domain on which we simulate our model represents only a small space around an individual islet. We note that this region is small compared to the entire pancreas, and further that the chemokine is free to diffuse out of our prescribed region. Typically, when solving PDEs, boundary conditions are applied at the edge of the domain over which the equation is being solved. However, as we are treating our region as only a small part of a much larger one, we feel it is more appropriate to not apply boundary conditions, leaving them open. The reaction term is given by local decay with point sources given by the locations of the beta cells. The chemokine equation now reads:

$D{\nabla }^{2}C=\lambda C-{\displaystyle \sum _{j\in \beta }{\nu }_j}(t)\delta (x-x_j),(4)$

where λ is the degradation rate of the chemokine and x_j are the centroids of the beta cells. To reflect the fact that dead beta cells will not secrete chemokine, we assume that the production of chemokine is dependent on the current viability of the cell ν_j(t). For a viable beta cell, ν = 1, whilst dead beta cells have ν = 0. This equation can itself be recast as:

$QC={\displaystyle \sum _{j\in \beta }{\nu }_j}(t)\delta (x-x_j),Q=(\lambda -D{\nabla }^2).(5)$

The Green's function for the operator Q at a given time is a Gaussian function and so we lump parameters together and choose a form for the chemokine signal given by:

$C(x,t)={\displaystyle \sum _{j\in \beta }\alpha }{\nu }_j(t)\exp \left(-{(x-x_j)}^2/(2{\sigma }_C^2)\right),(6)$

where α now represents the strength of the chemokine signal and σ_C represents the range over which it decays. In order to represent our assumption that the chemokine exists in forms in which it is bound to the membrane and forms in which it freely diffuses, we replace Equation (6) by:

$C(x,t)={\displaystyle \sum _{k=1}^2{\displaystyle \sum _{j\in \beta }\alpha }}{\nu }_j(t)C_0^k\exp \left(-{(x-x_j)}^2/(2{\sigma }_{C,k}^2)\right),(7)$

where σ_{C, 1} is small and reflects the membrane bound chemokine, whilst $C_0^k,k=1,2$ control the proportion of chemokine that is membrane bound vs. that which is freely diffusing. The immune cells respond to this gradient via:

$F_i^{\text{chemo}}=S_i\nabla C,$

where S_i is the sensitivity of the i'th cell to the chemokine gradient. We note that this represents the notion that immune cells tend to move up a potential gradient in which the locations of the beta cells are given by local peaks of this potential. The use of such potentials to model forces is an integral part of this modeling framework (Middleton et al., 2014).

3.5. Cell-Cell Interactions

We assume that T and B cells are attracted to one another. For simplicity's sake, we assume that this can be modeled via another Gaussian function. It is worth noting that this could represent attraction acting at a distance, or could represent a local effect acting to keep cells in contact with one another (or some combination of the two). We do not allow cells to pass through one another. We can achieve both of these through the use of another potential. Common choices for these potentials include the Lennard–Jones potential and the Morse potential (Middleton et al., 2014). We shall use the following potential:

$\begin{array}{l}{U}_i^{\text{cell-cell}}={\displaystyle \sum _{j\ne i,j\in \mathcal{T}\cup ℬ}{A}_{\text{a}}^{ij}\exp (-|x_j-x_i{|}^2/d_{\text{a}}^2)}\hfill \\ -A_{\text{r}}{\left(|x_j-x_i|-2r_{\text{imm}}\right)}^{6}H(|x_j-x_i|-2r_{\text{imm}}).\hfill & (8)\end{array}$

The first of these terms is the attraction between immune cells of different types—B cells seek out T cells to activate and T cells are also attracted to B cells. This is achieved by setting:

$A_{\text{a}}^{ij}=\{\begin{array}{ll}{A}_{\text{a}}^0\hfill & {\chi }_i\ne {\chi }_j\hfill \\ 0\hfill & {\chi }_i={\chi }_j\hfill \end{array},$

where χ_i is an indicator function taking value 1 if cell i is a T cell and 0 if it is B cell. The spatial scale of this attraction is set by d_a.

The second term accounts for the fact that cells cannot pass through one another. If the distance between cells is greater than the sum of their radii, that being 2r_imm, they exert no repulsive force on one another, achieved through the inclusion of a Heaviside function H. However, if this distance falls below this sum, a large repulsive force is exerted. This is known as a hard-core potential (Levine et al., 2000). Note that some authors choose this potential to be infinite at the point of contact, so that cells cannot occupy or share any part of physical space (Bruna and Chapman, 2012). However, this necessitates choosing a smaller time step and, since cells can remain in contact with one another for significant periods of time, we will not make such a choice. Note also that the parameter A_r is shared amongst all immune cells. Once again, the forces acting on the i'th cell can then be expressed as:

$F_i^{\text{cell-cell}}=-\nabla U_i^{\text{cell-cell}}.(9)$

Note that, since immune cells move around the domain, the potential “landscape,” together with its peaks and troughs, is constantly evolving.

3.6. Islet Interactions

For simplicity, and to remain consistent with our chosen geometry of the cells, the basement membrane encapsulating the islet is chosen to be a circle centred at the origin with radius R. Cells cannot pass through the membrane, so when it is present, immune cells remain on the inside or outside of the islet. This can be reflected through the use of another hard-core potential. Whilst we have not yet defined the dynamics for the membrane, we shall assume that there is a integrity threshold, h, below which the membrane does not prevent cells from passing through. Since we do not expect the membrane to have the same integrity across its whole length, we also need to account for the location at which the cell interacts with the islet. The relevant hard-core potential is then given by:

$\begin{array}{l}{U}_i^{\text{islet}}=A_{\text{islet}}H(m(\theta )-h)H(r_{\text{imm}}-|d_i-R|)\hfill \\ {\left(r_{\text{imm}}-|d_i-R|\right)}^6,\hfill \\ d_i=|x_i-x_{\text{islet}}|,\hfill & (10)\end{array}$

where m(θ) ∈ [0, 1] is the membrane integrity at angle $\theta \in {\mathbb{S}}^1$ around the membrane, $x_{\mathrm{\text{islet}}}\in {ℝ}^2$ is the location of the center of the islet and A_islet is a constant indicating the strength of repulsion. Once more, the force acting on the cells is then given by:

$F_i^{\text{islet}}=-\nabla U_i^{\text{islet}}.(11)$

Finally, we must define the forces acting upon the immune cells exerted by the beta cells. This is functionally the same as the repulsive force between cells (Equation 8) and is given by:

$\begin{array}{l}{U}_i^{\beta }=A_{\text{r}}{\displaystyle \sum _{j\in \beta }{\nu }_j}(t)H(|x_j-x_i|-(r_{\beta }+r_{\text{imm}}))\hfill \\ {\left(|x_j-x_i|-(r_{\beta }+r_{\text{imm}})\right)}^6,\hfill & (12)\end{array}$

where r_β is the radius of the beta cells. The time-dependence of ν_j here reflects that the beta cells may die, after which we no longer need to consider repulsive effects generated by them (assuming the dead cell bodies are cleared by macrophages), and so ν_j is set to 0 for that cell. The final force in Equation (2) is defined through

$F_i^{\beta }=-\nabla U_i^{\beta }.(13)$

3.7. Immune Cell Lifespan

The timescale over which insulitis takes place is long compared to the average lifespan of an immune cell. Thus, we need to incorporate immune cell death into our model. We assume that T and B cells have a lifespan of L_T, and L_B respectively. Each cell then has a counter χ_i, which is incremented by setting χ_i ↦ χ_i + 1 at each time step. Activated T cells are expected to have shorter lifespans than unactivated T cells (Green et al., 2003), so for those cells, we instead update the counter via χ_i ↦ χ_i + Δχ where Δχ > 1.

When χ_i exceeds L_T (or L_B) for a given T (B) cell, it is considered to be dead and is removed. We keep the number of T and B cells constant throughout the simulation by assuming that each dead cell is replaced by a newly arriving one. We assume that the vasculature is sufficiently dense that immune cells can enter at any point in the extra-islet space. As such, the location of the new cell is drawn randomly from a uniform distribution over the extra-islet space, achieved through the use of polar coordinates. The cell's counter is reset to 0 and, if the cell is a T cell, it is chosen to be in the unactivated state.

3.8. Membrane Dynamics

In our simplified geometry, the encapsulating membrane is a represented by a circle, which can be parameterized by a single variable $\theta \in {\mathbb{S}}^1$. The equation governing the evolution of the membrane viability, m, is given by:

$\begin{array}{l}\frac{\text{d}m(\theta )}{\text{d}t}={\alpha }_m(1-m(\theta ))\hfill \\ -{\lambda }_m{\displaystyle \sum _{i\in \mathcal{T}}\exp }\left(-(|x_i-x_m(\theta )|-r_{\text{imm}})/d_m\right),\hfill & (14)\end{array}$

where α_m and λ_m are respectively the repair and degradation rates of the membrane and x_m(θ) is the location in real space of the membrane at position θ. The term in the sum represents the fact that we expect T cells to break down the membrane and this can only occur when T cells are close enough. The range over which the T cells can degrade the membrane is set by d_m.

We note that as beta cells are destroyed, the morphology of the affected islets may change. In particular, some islets in T1D present with small size and irregular outline (Gepts, 1965). The membrane in our study, where it exists, simply follows a circle with fixed radius, and thus does not account for these changes. However, we note that not all degranulated islets have small size and irregular outline (Gepts, 1965). Moreover, the morphology of islets, irrespective of insulitis, is highly variable. It is thus difficult to infer how the morphology of a given islet varies over the course of the insulitic process. Nevertheless, it seems likely that changes to islet morphology will occur, and we discuss approaches to incorporate this into our model in Section 5.1.1.

3.9. Activation and Apoptosis

When T cells and B cells are in contact with one another, we assume that B cells can activate non-activated T cells and that this process has a characteristic time course. We thus describe the activation level of the i'th T cell, a_i, via:

$\frac{\text{d}{a}_i}{\text{d}t}=\{\begin{array}{ll}\begin{array}{l}\sum _{j\in ℬ}H(|x_j-x_i-2r_{\text{imm}}+0.1|)\hfill \\ -{\lambda }_a{a}_i\hfill \end{array}\hfill & a_i<1\hfill \\ 0\hfill & a_i\ge 1,\hfill \end{array},i\in \mathcal{T},(15)$

where λ_a is the decay rate of the activation signal. The constant 0.1 is included inside the Heaviside function to account for the fact that we are not explicitly modeling cell-cell contacts, and the use of hard-core potentials tends to make cells move apart quickly when they are close. Our specific choice for this constant is selected such that cells are allowed to remain in contact with one another. We note that if this value is selected to be too large, cells can essentially occupy the same location, whilst values that are too small will cause cells to “bounce” off one another. We wish to avoid both of these behaviors, and over a range of choices, we found that a value of 0.1 satisfied these criteria. Once the activation signal reaches 1, the cell is activated and remains so until that cell dies, (or until the end of the simulation if that occurs first), reflected by setting the RHS to zero upon the activation reaching 1. Activated cells have an increased sensitivity to the chemokine signal and enhanced killing rate. This is captured in the model by increasing the sensitivity and killing parameters:

$S_i\mapsto S_i+\Delta S,{\kappa }_i\mapsto {\kappa }_i+\Delta \kappa ,$

where ΔS and Δκ are positive constants.

The dynamics for the apoptosis of beta cells b_i follows a similar prescription:

$\frac{\text{d}{b}_i}{\text{d}t}=\{\begin{array}{ll}\begin{array}{l}\sum _{j\in \mathcal{T}}{\kappa }_{j}H(|x_j-x_i-2r_{\text{imm}}+0.1|)\hfill \\ -{\lambda }_b{b}_i\hfill \end{array}\hfill & b_i<1\hfill \\ 0\hfill & b_i\ge 1,\hfill \end{array},i\in \beta ,$

where λ_b is the decay rate of the apoptotic signal, κ_j is the killing rate of the j'th T cell, and we include the constant, 0.1, for the same reasons as for (Equation 15). If b_i exceeds 1 for a given cell i, that cell is assumed to have been killed and is removed from the simulation. This is achieved by setting its viability, ν_i, to be zero.

3.10. Implementation

Since the system we are solving is a Langevin equation, we must use an appropriate numerical method to deal with the stochasticity. For computational efficiency, we shall use the forward Euler–Maruyama scheme. This method only has strong order dt^1/2 but requires the fewest function evaluations of any of the solvers for stochastic differential equations (SDEs). The one exception to this is the equation governing the membrane dynamics where the explicit Euler scheme is unstable and causes solutions to blow up. For this equation, we use an implicit, backward Euler scheme. This keeps the total order accuracy of solutions consistent with the other equations, but provides the required stability without incurring additional computational cost.

Note that the membrane (Equation 14) treats distinct points along the membrane as being distinct from one another. This is the only system that is grid-based – all of the other equations are grid-free. In order to solve this system, we first discretize the membrane into N_θ points, θ_i = − π + 2π(i − 1)/N_θ, i = 1, N, and subsequently solve (Equation 14) at each of these points. To find the value of the m at a point not on the discrete grid, we use band-limited interpolation for periodic signals as described in Schanze (1995).

3.11. Initial Conditions and Parameters

The immune cells are initially located at random, non-overlapping positions within annular domain with inner radius R and outer radius 1mm. The islet is positioned at the center of the domain, which for simplicity, we set as our origin. To distribute the beta cells, N_β cells are arranged in concentric circles within the islet. Next, for T = 1000 timesteps, they are allowed to freely evolve governed by:

$\frac{\text{d}{x}_i}{\text{d}t}=F_i^{\text{islet}}+F_i^{\beta -\beta }+\sigma {\xi }_i(t),i\in \beta ,(16)$

where x_i is the center of the i'th beta cell, $F_i^{\mathrm{\text{islet}}}$ is the same as Equation (13) operating on the beta cells rather than the immune cells with:

$F_i^{\beta -\beta }=-\nabla U_i^{\beta -\beta },(17)$

where

$U_i^{\beta -\beta }={\displaystyle \sum _{j\ne i,j\in \beta }H}(|x_j-x_i|-2r_{\beta }){(|x_j-x_i|-2r_{\beta })}^6,(18)$

where $x_i\in {ℝ}^2$ now represents the location of the beta cell. After the n steps, the beta cells are fixed in location for the remainder of the simulation. Note that during this process, the immune cells are fixed in position. We distribute cells in this manner as it (a) generates a more realistic islet geometry compared with spacing them evenly within the islet and (b) it overcomes the known computational problems with randomly distributing non-overlapping disks in a confined region (Song et al., 2008). In Figure 1, we show the initial and final configuration of the beta cells following this approach. All simulations and analyses were performed in Matlab.

FIGURE 1

Figure 1. (A) Initial configuration of beta cells within an islet, distributed in concentric circles. (B) Final configuration of the beta cells evolved for T = 1000 timesteps, according to Equations (16–18).

Other initial conditions are given by m = 1 for all θ, a_i = 0 for all $i\in {T}$, and b_i = 0, ν_i = 1, ∀i ∈ β.

The sensitivities of the T and B cells to the chemokine gradient are drawn from distinct Gaussian distribution with means S_T and S_B and variances σ_T and σ_B respectively. Parameters governing the geometry of our domain, such as immune cell and beta cell sizes were chosen to be matched to available data. In particular, beta cell diameters were taken from Saisho et al. (2013), immune cell diameters were based on values found in Wang et al. (2012), numbers of cells in each immune cell population were taken from Willcox et al. (2009) and islet composition (in terms of beta cells) were based on those found for humans in Kim et al. (2009). Data regarding the in vivo lifespan of specific effector lymphocytes in humans are rare and these may also be heterogeneous even within the same population. However, effector lymphocytes are unlikely to have lifespans ranging beyond several weeks (Sprent, 1993). To address the uncertainty in lifespans, we performed simulations across a range of values, from 3 to 7 weeks and compared results.

Other parameters in our model were tuned to provide the required behaviour, though we comment that much of the general behavior observed in our simulations is robust to parameter variations, suggesting that the model outcomes are robust. We note that certain processes, such as the degradation and repair of the peri-islet basement membrane, are phenomenological in nature, and as such, it is difficult to ascribe to them meaningful parameter values. Whilst the islet membrane is a physical component of the real system, we represents its viability by a scalar in the range 0–1, and so parameter values should be interpreted with respect to this scaling. Similar comments hold for the chemokine signal. Though a chemoattractant gradient is thought to exist in the peri-islet space, it is not clear what its composition, and subsequent properties, might be. As such, we have chosen a prototypical form to represent our chemokine signal, and investigated how changes to the strength of this signal impacts upon the resulting dynamics. All parameter values are summarized in Table 1 and referenced in text where they are altered for specific numerical experiments.

TABLE 1

Table 1. Table of parameter values and meanings for the agent-based model.

4. Results

In Figure 2, we show the initial state of a prototypical simulation, and the state of the system at T = 4000 days. Here, we can see that the islet exhibits significant beta cell loss. We also see that immune cells congregate at the islet membrane and that a small number infiltrate the islet itself. Across all parameter sets used in our experiments, this behavior is preserved. In the Supplementary Material, we provide videos showing the full time-course of typical simulations with base parameters as indicated in Table 1, lymphocyte lifespan of 21 days and other parameters selected to reflect a variety of conditions, as discussed in the forthcoming sections. These simulations may be compared with the averaged results shown in Figures 3–7.

FIGURE 2

Figure 2. (A) Initial state of system. (B) State of the system after T = 4000 days. A video of the full time-course of this simulation can be found in the Supplementary Material. In both panels, the blue circles represent beta cells, the red circles represent T cells, the green circles represent B cells, the black circle represents the peri-islet basement membrane and the black gradient represents a chemokine signal. Parameters are as in Table 1.

FIGURE 3

Figure 3. Effect of variation of B cell number on beta cell mass loss in the agent-based model. Each graph shows the mean ± SEM for 100 simulations. Across all immune cell lifespans, we find that increasing the number of B cells significantly increases the destruction rate of the beta cells. Parameters are as in Table 1. Videos of typical simulations for the case with immune cell lifespans of 21 days may be found in Videos S1, S2.

For all parameter sets, we simulate 100 realizations and average at specific time points to develop a time-course. During the simulations, we keep track of the number of immune cells within the islet and the number of viable beta cells. We take the total area encompassing the viable beta cells as being reflective of the beta cell mass. We shall now report both the time-courses of the remaining beta cell mass as a percentage of the initial mass (calculated as the ratio of the total mass of viable beta cells at time t and the initial mass of beta cells, averaged over the 100 realizations), and the number of immune cells in the islet for a number of conditions, along with the standard error of the mean (SEM) at each time point. In each case, we shall use the base set of parameters as indicated in Table 1, changing specific values where indicated. In all cases, to assess how the lifespan of the immune cells affects the resulting dynamics, we vary this between 21, 28, and 56 days as indicated in the figures.

4.1. Varying Number of B Cells

We begin by examining a key condition in this study, namely how the number of B cells in the peri-islet space impacts the rate of destruction of beta cells. To address this, we compare a condition with few B cells, setting N_B = 5, and one with a higher number of B cells, with N_B = 30. These results are summarized in Figure 3.

We observe that when the immune cell lifespan is short, neither the low nor the high B cell condition exhibit significant beta cell destruction. As the lifespan is increased, we observe a pronounced increase in this rate. For an immune cell lifespan of 56 days, we see that the high B cell condition reaches an equilibrium value at around 3000 days. Note that the equilibrium value is below the critical threshold of ~20% of remaining beta cell mass (Cnop et al., 2005) (however, it should be noted that this threshold is a whole-body threshold and may not be valid for an individual islet).

4.2. Varying Repair Rate of Membrane

We have assumed in our model that the peri-islet basement membrane is capable of being repaired at a rate α_m. If we decrease α_m, the rate of degradation of the membrane may be such that holes caused by infiltrating T cells are not repaired (at least over the timescale of our simulation). In the extreme case, the membrane is incapable of repair, in which case we set α_m = 0. In Figure 4, we display the time-course of a simulation in such a regime.

FIGURE 4

Figure 4. Time-course of beta cell mass comparing conditions with and without peri-islet basement repair. Each graph shows the mean ± SEM for 100 simulations. Data for simulations with 5 B cells are displayed in the left column, whilst those for simulations with 30 B cells are depicted in the right column. In the case without membrane repair, we observe both an increase in the rate of beta cell destruction and pronounced differences in the final beta cell mass across all immune cell lifespans and B cell number. Parameters are as in Table 1 with α_m = 0. Videos of typical simulations for the case with immune cell lifespans of 21 days may be found in Videos S3, S4.

It is clear from these graphs that removing the repair mechanism of the membrane dramatically increase the susceptibility of the beta cells to cytotoxicity. In the most extreme case, in the situation with few B cells and a short immune cell lifespan, we now see a huge difference in both the rate of beta cell destruction and the final cell mass at the end of our simulation. In the case with no membrane repair, the final beta cell mass of almost all simulations is below ~20% of the initial mass. These results highlight the importance of the peri-islet basement membrane to the progression of insulitis, and also suggest that ongoing repair of this membrane is critical to slowing the infiltration of the immune cells.

We also observe, in contrast to the results in Figure 3, that the immune cell lifespan now has a more muted affect on the resulting dynamics. Whilst there is an increase in the rate of beta cell destruction with increasing lifespan, the more significant contribution appears to be to lower the equilibrium value of remaining beta cell mass. This likely arises due to the fact that, as the beta cell population decreases, so too does the chemokine signal to which the immune cells respond. This means that the cells (on average) move into the islet at a slower rate, since their paths becomes more dominated by noise. If the cell lifespan is too short, the immune cells may die before entering the islet if there are insufficient beta cells generating a chemokine gradient.

4.3. Increasing the Membrane Degradation Rate

Instead of reducing the basement membrane repair rate, we now increase the efficacy of the immune cells of degrading the membrane, setting λ_m = 0.2 and display the results in Figure 5. For comparison with the condition in which the encapsulating membrane is not repaired, we additionally plot results from this case in the same figure, as indicated.

FIGURE 5

Figure 5. Time-course of beta cell mass comparing with higher membrane degradation rates. Each graph shows the mean ± SEM for 100 simulations. The left column shows results for simulations with 5 B cells whilst the right column shows those with 30 B cells. Data for the condition with the base value for membrane degradation but with α_m = 0, as reported in Figure 4 are also shown. An increase in the beta cell destruction rate, as well as lower final beta cell masses are observed across all immune cell lifespans and number of B cells. Note, however, that whilst the initial rate of beta cell destruction for the condition with high membrane degradation is higher than that for the case with no membrane repair, the final beta cell mass is higher for the former compared with the latter. We refer the reader to the text in Section 4.3 for further interpretation of this observation. Parameters are as in Table 1 with λ_m = 0.2 for the high degradation condition and α_m = 0 for the no membrane repair condition. Videos of typical simulations for the case with immune cell lifespans of 21 days may be found in Videos S5, S6.

As expected, increasing the rate of membrane degradation increases the rate of beta cell destruction and subsequently decreases the final remaining beta cell mass. Compared with the condition with no membrane repair, we see that the initial rate of beta cell destruction is greater in the case with high membrane degradation. However, as the insulitis process continues, this trend becomes reversed. Remarkably, across all cases varying the number of B cell and the immune cell lifespan, the time at which this occurs appears to be approximately the same. We also see that the final beta cell mass is lower for the condition with no repair than for that with high membrane degradation.

At the beginning of simulations, the ability of the T cells to degrade the basement membrane at a faster rate allows them to infiltrate the islet and subsequently kill the beta cells more quickly. Whilst the initial degradation of the membrane in the no repair condition is slower, the fact that “holes” in the encapsulating membrane caused by T cells are not repaired means that, at later times, immune cells can simply pass through these holes.

When the membrane can be repaired, the rate of immune cell infiltration is governed by the balance between how quickly it can be degraded and how quickly it can be repaired. As the chemokine signal across the peri-islet space becomes weaker when beta cells die, the average rate that T cells reach the basement membrane falls, since the strength of the signal is a contributory factor in determinining the transit time for a T cell toward the islet. As there are fewer T cells aggregrating at the islet membrane at any given time, the degradation of the membrane becomes slower, and hence the beta cells are destroyed at a slower rate as time progresses. At later times the arrival rate of T cells may be so slow that the membrane is repaired at the same rate as it is degraded, on average, and so islet infiltration by the immune cells will be halted.

These observations account for the slower rate of beta cell destruction and higher remaining beta cell mass for the high degradation condition compared to that with no membrane repair. This further highlights the importance of continual basement membrane repair to the slowing of insulitis.

4.4. Increasing the Killing Efficiency of the T Cells

We now increase the killing rate of the T cells. To reflect the fact that we assume that activated T cells are responsible for most of the beta cell destruction, we achieve this by doubling the killing rate boost, setting Δκ = 1.8. The results for these experiments are displayed in Figure 6.

FIGURE 6

Figure 6. Time-course of beta cell mass comparing with higher activated killing rates. Each graph shows the mean ± SEM for 100 simulations. The left column shows results for simulations with 5 B cells whilst the right column shows those with 30 B cells. Although an increase in the rate of beta cell destruction is observed, this is relatively small in all cases. Parameters are as in Table 1 with Δκ = 1.8.

We observe that increasing the killing rate of the immune cells does not significantly affect the death rate of beta cells, either quantitatively or qualitatively. In all cases, we observe an increase in the beta cell mass destruction rate, but this is fairly minimal. The limited impact of increased killing rate can be understood by noting that, whilst the killing rate is an important factor in the overall insulitis process, this is dominated by transit time to and within the islet, as well as activation by B cells. Once inside the islet, T cells can interact with and kill beta cells but the total number of beta cells an individual T cell can destroy will be influenced more strongly by its remaining lifespan than on its killing rate. This highlights the need to incorporate the spatial aspects of the immune response and that high T cell killing rates do not necessarily significantly affect the progression of insulitis.

We comment that these observations at first seem at odds with previous results from the NOD mouse where the avidity of a T cell population is well-correlated with its pathogenic potential, and the progression of T1D occurs through avidity maturation (Amrani et al., 2000). Theoretical studies have supported the viewpoint that higher avidity CD8+ T cell populations lead to increased rates of beta cell destruction (Khadra et al., 2011) and that low avidity populations may have a protective effect against T1D development (Khadra et al., 2009). However, recent results in the NOD model have suggested that the highest avidity CD8+ cells do not escape thymocyte negative selection, and accordingly, the cells responsible for beta cell destruction in the murine pancreas may be regarded as those with “intermediate avidity” (Han et al., 2016). Furthermore, the paucity of data on the specificity of lymphocytes involved in human insulitis makes it difficult to propose quantitative statements about the effective killing potential of individual T cells. This is particularly difficult since the majority of the results regarding the avidity of human lymphocytes come from peripheral blood measurements and are not necessarily indicative of avidities within the pancreas (Reijonen et al., 2004; Standifer et al., 2009). Finally, the lifespan of a T cell that has successfully infiltrated an islet should still have a significant impact on the number of beta cells it can destroy, assuming its avidity is sufficiently high.

4.5. Increasing the Rate of Chemokine Production

To reflect increases in the rate of chemokine production, we now double ${\alpha }_C^1$ and ${\alpha }_C^2$ simultaneously, setting ${\alpha }_C^1=1.0$ and ${\alpha }_C^2=20.0$. Note that, since we are assuming a quasi-steady state approximation for the chemokine, this simply adjusts the profile of the chemokine signal across our domain. The results for these simulations are depicted in Figure 7.

FIGURE 7

Figure 7. Time-course of beta cell mass comparing conditions with the base and increased strength of chemokine signal. Each graph shows the mean ± SEM for 100 simulations. The left column shows results for simulations with 5 B cells whilst the right column shows those with 30 B cells. In the high chemokine signal parameter regime, we observe both an increase in the rate of beta cell destruction and pronounced differences in the final beta cell mass across all immune cell lifespans and B cell number. These differences are greater than those displayed for higher T cell killing rates in this figure. Parameters are as in Table 1 with ${\alpha }_C^1=1.0$ and ${\alpha }_C^2=20.0$. Videos of typical simulations for the case with immune cell lifespans of 21 days may be found in Videos S7, S8.

As for the situation with increased killing rates, we observe similar qualitative behavior to the case with lower chemokine production rates. However, in contrast to the former case, the quantitative differences between the beta cell death rate and final cell mass are much greater. This further highlights the contribution of cell movement in the overall inflammatory response.

4.6. Combining High Chemokine Signals and Reduced Membrane Repair

In our final numerical experiment, whose results are shown in Figure 8, we simultaneously increase the strength of the chemokine signal, setting ${\alpha }_C^1=1.0,{\alpha }_C^2=20.0$, as in the previous subsection, and remove the membrane repair, setting α_m = 0. As expected, with these combined alterations, we now observe significantly high rates of beta cell destruction for all cases considered. Moreover, the final beta cell mass at the end of our simulations is below 20% of the original cell mass. For the cases with the longest lived immune cells, those having a lifespan of 56 days, we now see that essentially all of the beta cells have been destroyed within 1000 days.

FIGURE 8

Figure 8. Time-course of beta cell mass with high chemokine signals without basement membrane repair. Each graph shows the mean ± SEM for 100 simulations. The left column shows results for simulations with 5 B cells whilst the right column shows those with 30 B cells. In all cases, the rate of beta cell loss is high and the remaining cell mass at the end of the experiments is low. Parameters are as in Figure 7 with α_m = 0.

4.7. Immune Cell Invasion Profiles

We now consider how the invasion of immune cells is dependent on the remaining beta cell mass and report the dependence of the number of immune cells within the islet (i.e., within the basement membrane) on the remaining beta cell mass. These results are shown in Figure 9 for the case with 5 B cells and in Figure 10 for the case with 30 B cells. The results for the case with high chemokine signals with no basement membrane repair, and for the case with increased membrane degradation are omitted, but are qualitatively similar to the case with no repair alone. In both figures, we only show results when the immune cell lifespan is 56 days.

FIGURE 9

Figure 9. Profiles of immune cell invasion dependent on the remaining beta cell mass in cases with 5 B cells. Each graph shows the mean ± SEM for 100 simulations. We report the number of T and B cells within the islet in a number of conditions. The top-left panel shows results using the parameters in Table 1, the bottom-left uses these values with Δκ = 1.8, in the top-right, we set ${\alpha }_C^1=1.0$ and ${\alpha }_C^2=20.0$, and the bottom-right has α_m = 0.

FIGURE 10

Figure 10. Profiles of immune cell invasion dependent on the remaining beta cell mass in cases with 30 cells. Each graph shows the mean ± SEM for 100 simulations. We report the number of T and B cells in the same conditions as displayed in Figure 9.

For situations with 5 B cells, we observe two qualitative behaviors amongst the range of cases considered. Using the default parameters and those with high killing rates, we observe no infiltration of the islet by the B cells. The immune cell population undergoes an initial increase to a maximum value, followed by an approximately linear decay with the remaining beta cell mass. The behavior with high chemokine signals or no membrane repair does show islet invasion by the B cells. In these cases, the T cell population now both increases and decreases more gradually as a function of the beta cell mass, and attains a higher maximum value than in either of the two previous cases.

If we now increase the number of B cells to 30, we observe, in all cases, invasion of the islets by B cells and a qualitatively similar profile of T cell infiltration. With the default parameters, or those with high killing rates, the maximum value of invading T cells occurs when approximately 50% of the beta cell mass remains. In contrast to the T cell population, the number of invading B cells varies little with the beta cell mass. In the case with no membrane repair, this is no longer true, and we now find that the B cell profile follows the T cell profile and furthermore, that the maximum value of immune cell infiltration occurs at low values of the remaining beta cell mass.

4.8. Trial-to-Trial Variability and the Role of Space

For each of the scenarios listed, we report the standard error of the mean for all recorded values at each time point. In general, we observe that variability across our simulations for a given scenario is low, suggesting that the model behavior is robust and that noise does not significantly affect the intrinsic model dynamics. In particular, scenarios with either very high rates or very low rates of beta cell destruction exhibit very little variability. Variation is highest in scenarios that lead to intermediate rates of beta cell death, suggesting that the role of noise in such situations is much more significant in these conditions.

The results reported in Figures 3–10 reflect spatial averages that are taken over the peri-islet space. As such, they do not reflect the spatial component of the dynamics produced by the agent-based model. These effects are, in general, quite subtle. In particular, since we assume that activation of CD8+ cells by CD20+ cells requires them to be co-localized, the spatial arrangement of the immune cells will play an important role in determining the overall rate of beta cell destruction. The geometry in this initial version of the model has been chosen to be simplistic in nature. In future versions, the geometry will be adjusted to better match that of real islets as discussed in Section 5.1.1. One could imagine scenarios in which the routes taken by immune cells when migrating toward the islet could influence the average number of interactions with one another, and further, that these paths would be influenced by prominent vasculature near the islet.

The spatial configuration of the immune cells also plays a significant role in determining how quickly the basement membrane is degraded. Where many CD8+ cells accumulate at specific points along the periphery of the islet, we would expect high rates of membrane degradation. Conversely, if the cells are evenly distributed along the membrane, average degradation rates are likely to be low. Our results suggest that the dynamical processes governing the repair and destruction of the basement membrane play a key role in determining the overall rate of beta cell destruction and so it is clear that the spatial distribution of immune cells along this membrane will also be critical.

These spatial considerations are the key determinant for the degree of variability observed between trials, as reported in each of the figures. In the Supplementary Material, we provide animations of typical simulations in a subset of the scenarios listed to demonstrate the subtle role that these effects have on the insulitis process.

5. Discussion

In this manuscript, we have constructed a spatial, agent-based model of immune cell invasion of a prototypical islet of Langerhans. The agents in the model are immune and beta cells, and each obeys rules that broadly match the general behavior of these cell types. The model attempts to mimic a real islet, allowing for in silico experiments to be performed through model simulation. The agent-based framework allows both the spatial aspects of insulitis to be investigated, and avoids issues of applying ODE-based approaches to a system that has low cell numbers. It is intended that this will lay down the framework for more extensive spatio-temporal modeling of insulitis, with an ultimate aim to understand its core mechanisms and to devise strategies to slow or halt its progress.

Overall, the observed phenomena in the model are well matched to real data, which have been collected and analyzed from human patient samples (Willcox et al., 2009; Morgan et al., 2014), and are reported in Figure 11. We note that data for the pauci-immune phenotype have recently appeared in Leete et al. (2016), but since our model is primarily focussed on the hyper-immune phenotype, we highlight these results here. Results from the agent-based model have the same qualitative behavior as displayed by these data. This in turn provides evidence that the mechanisms included within our model are sufficient to capture qualitatively the key features of the inflammatory response. In particular, we note that the case with no membrane repair appears to provide the best qualitative match to the existing data, in which we see a rise in the number of infiltrating immune cells during disease progression, with a peak at around 10–20% of remaining beta cell mass. We note that scenarios with low overall rates of beta cell death are associated with low levels of B cell infiltration. This further supports the viewpoint that the CD20+ cells play an important role in determining the rate of beta cell destruction. By adding and removing mechanisms in a systematic way, we can probe what effect their inclusion or exclusion will have on the resulting dynamics and subsequent prognosis. It is important to note that, since we investigate the disease at the level of the islets, we can make predictions about how modifying certain cell properties will affect the inflammation of the islets, rather than considering whole body responses.

FIGURE 11

Figure 11. Data showing the immune cell population against the percentage of the islet which is positive for insulin, as calculated via morphometry. This percentage can be thought of as being representative of the residual beta cell mass. Each point is the mean ± SEM for the number of relevant immune cells detected in each islet section, averaged over 279 islets across 29 patients with T1D. The number of islets analyzed in each insulin-positive category is given in parenthesis. Permission to reprint this figure from Willcox et al. (2009) was granted by John Wiley and Sons (Licence number: 3971831129485; www.interscience.wiley.com).

Accordingly, in the present study, we have performed several parameter studies to examine how the relative contributions of core processes affects insulitis. We have demonstrated that, under the assumptions used in constructing the model, an increased number of B cells in the peri-islet space gives rise to a faster destruction of the beta cells. This effect can be amplified by factoring in manipulations of other processes, such as the killing rate of the beta cells by T lymphocytes and the repair rate of the islet's basement membrane. Simply increasing the killing rate of the T cells did not have a significant impact of the rate of the beta cell destruction. However, increasing the strength of the chemokine signal to which the immune cells respond did show a pronounced increase in the death rate. These facts highlight the importance of the chemotactic process by which the immune cells enter the islet, find and kill the beta cells.

In scenarios in which the islet basement membrane is unable to repair itself, or is degraded more rapidly, we also see significant increases in the beta cell destruction. Moreover, we also find that the final beta cell mass at the end of our simulations is substantially depleted. This highlights the importance of this membrane as a barrier to the infiltrating immune cells. In particular, when chemokine signals are sufficiently high, if the membrane is incapable of being repaired, all of the beta cells may ultimately be killed, whereas in other situations, we find residual beta cells even after long simulation times (not shown). The importance of the basement membrane during insulitis in our model is in line with a recent study of this membrane in samples from human patients with T1D (Bogdani, 2016).

Upon examining the dependence of the immune cell invasion on the remaining beta cell mass, we find in most cases the profile of B cell infiltration is relatively flat. In cases with low B cell numbers, low chemokine signals and with a membrane that is repaired, there is essentially no B cell invasion. In situations with higher numbers of B cells, the beta cell mass at which the T cell invasion is maximal is shifted to lower values. In the case with no membrane repair, this maximal point is shifted to very low values of the remaining cell mass. Moreover, in this instance, we also see that the B cell invasion profile follows that of the T cells. This condition is most consistent with the human data reported in Willcox et al. (2009), suggesting that, for our model, low membrane repair rates are an important factor in the insulitis process.

The rise and subsequent decrease in T cell populations observed in our simulations is in accordance with experimental results of T1D specific autoimmune responses in the NOD model (Trudeau et al., 2003). These dynamics have been described as “waves” and have been explained mathematically by transitions through Hopf and homoclinic bifurcations in Mahaffy and Edelstein-Keshet (2007) and by transient bistability of the autoimmune state in Jaberi-Douraki et al. (2014b). Interestingly, in the latter of these two theoretical studies, inclusion of multiple clones of T cells with different avidities can give rise to multiple waves in the T cell dynamics, which can significantly impact the rate of overall beta cell destruction. This highlights potential avenues for intervention through the promotion of lower avidity T cell populations.

It should be noted that data on the specificities of T cells within human insulitic islets are rare, but that cells with varying antigen specificities may be present (Coppieters et al., 2012; Babon et al., 2016). Thus, a potential development of our model would be to include multiple clonal populations to explore these dynamics in the spatial context, particularly given the potential insights that may be applicable from the theoretical results in Mahaffy and Edelstein-Keshet (2007) and Jaberi-Douraki et al. (2014b).

Agent-based and lattice-based approaches are becoming increasingly popular in biological modeling. ODE based approaches, which have long been a linchpin of the mathematical biology community become ill fit for purpose when wishing to describe processes in which the number of agents is small. Considering stochastic variants of ODEs and partial differential equations (PDEs) models is one potential approach to account for the variability that arises from low cell numbers, but this too is only a good approximation when the number of cells is sufficiently high. Moreover, these approaches do not typically account for trial-to-trial variability that may be important for describing individual results. Note that, at various levels of description, systems of differential equations can still form part of agent-based models.

One of the attractive properties of ODE and PDE modeling when compared with agent-based simulations is the computational cost associated with the latter, whereas density (continuum) based methods typically are low dimensional in nature. Indeed, a significant barrier to the widespread use of agent- and lattice-based models in the past has been the requirement of significant computational power. In recent years, the technological advances in computing have started to erode this barrier. In addition, domain decomposition methods (Chen et al., 2007; Tapia and D'souza, 2011) and the ability to parallelise computations on large scales on commercially available graphics processing units (Chen et al., 2007; Tapia and D'souza, 2011; Harvey et al., 2015) has made agent-based simulations achievable universally.

Here, we have focussed on a specific type of agent-based model. This is based on a discrete particle formulation, where each cell is parameterized by a finite area and location. Broadly, agent-based models can be divided into lattice-based and lattice-free methods (Plank and Simpson, 2012). The discrete particle approach is lattice-free, meaning that cells are free to move anywhere in the domain. One limitation of this approach is that the cells are restricted to be described as circles; that is, their shape cannot be deformed. Immune cell deformation may play an important role when considering how cells pass through degraded parts of the basement membrane. In the current model formulation, cells must “wait” until the gaps in the membrane are sufficiently large before they can pass through.

An alternative to the discrete particle formulation is the lattice-based cellular Potts model. This modeling framework has been growing in popularity and has been applied to cell migration in a number of contexts, including vasculogenesis (Merks et al., 2006; Daub and Merks, 2013), morphogenesis (Chen et al., 2007; Marée et al., 2007) and wound healing (Scianna, 2015). Although the cellular Potts model is lattice-based, the cells within the framework can occupy many lattice sites and thus their shape, as well as their location, evolves over time. On a fine enough lattice, such models can look remarkably similar to the system they are approximating. However, with the increase in spatial resolution of the lattice comes an associated computational cost. It should be noted that there are several actively maintained packages for agent-based modeling, such as CHASTE (Mirams et al., 2013) and CompuCell3D (Swat et al., 2012).

One approach that combines the flexibility of lattice-free approaches together whilst allowing for cells to change shape is to treat points on the cell membrane as particles and resolve forces acting upon them, under the assumption that the cell remains bounded by its membrane (Elliott et al., 2012). However, the mathematical machinery and computing power required to embed such a description into our domain are prohibitive for our purposes. A “halfway” house could, in which the cells are treating as deformable ellipsoids, is another possible way of relaxing our assumption that cells are perfectly circular (Palsson and Othmer, 2000).

5.1. Future Directions

The modeling framework described within this manuscript is not intended to be a fully comprehensive endpoint for the study of the spatio-temporal dynamics associated with insulitis. Instead, it is intended to be the first step to building a general virtual environment in which to simulate the invasion of islets during T1D. We will now highlight some potential areas for model development.

5.1.1. Geometry

We have treated perhaps the simplest possible geometry that is reflective of the islet. The model is planar in nature, whereas real islets and immune cells are three dimensional structures. We have ignored the mechanical effects on the basement membrane: as beta cells are removed from the islet, the morphology of the islet itself may change (Brereton et al., 2014). There is currently no structure imposed upon the extra-islet space. Accurate modeling, for example, in which the vasculature is taken into account may constrain locations at which immune cells are allowed to enter the domain. Through the use of imaging techniques, it should be possible to build a more faithful representation of both the intra- and extra-islet space and include these in the model.

5.1.2. Different Cell Types

The model could be developed to include more cell types, such as the macrophages that clear apoptotic beta cells, or different endocrine and exocrine cells which are not targeted by immune cells, but may slow their movement. There may be other cell-cell interactions that may be important during insulitis, for example the role of CD4+ cells in the activation of CD8+ cells (Castiglioni et al., 2005). Even restricting the model to the two immune cell types described in this manuscript, we could consider different clonal populations of the T cells (Khadra et al., 2011). One hypothesized strategy to slow the progression of insulitis is to promote the replication of a low avidity clone of T cells, that will out compete the higher avidity clones for space and resources, leaving only relatively passive T cells (Amrani et al., 2000; Khadra et al., 2009; Bluestone et al., 2015). Such a mechanism could be tested in our modeling framework.

5.1.3. Cell Proliferation

Currently, none of the cell types included in the model undergo cell division. In general, post-natal beta cell replication in humans is thought to occur at low rates, and we have disregarded it here. Similarly, immune cell replication in and around the islet is also assumed to be negligible. However, there is some evidence to suggest both that beta cell proliferation is enhanced during insulitis, at least during its early phases (Willcox et al., 2010; Dirice et al., 2014) and that immune cell replication also takes place within the islet itself (Graham et al., 2012) in the NOD model, though these findings have not been corroborated in human tissue. In particular, data collected from ex vivo human samples suggested that very few lymphocytes in the peri-islet space were positive for proliferative markers (Willcox et al., 2010). In order for the hypothesis of promotion of low avidity clones to be tested, the model would have to include cell division, since the action of long-lived memory cells is an integral part of that theory, however this may require coupling of the agent-based model with one representing the lymphatic system.

5.1.4. Multiple Islets

One of the striking features of insulitis is its heterogeneity, even within an individual patient. Islets near to one another can display different inflammatory profiles—insulitis in one islet thus does not imply insulitis in nearby islets. By extending our domain to include multiple islets, we can use the modeling framework to investigate under what scenarios this can occur. In particular, we can hope to identify what cellular interactions, on a broad scale, must occur for such heterogeneity in disease progression to be observed. One simple way to incorporate the influence of other islets in the current framework would be to include additional point sources of chemokine outside of the working domain. This would have subtle effects on the dynamics of cell movement, particularly when the number of beta cells becomes low.

5.1.5. Membrane Degradation

In the current model prescription, the membrane integrity is assumed to be compromised in a simple way by the action of nearby immune cells. This process is solely dependent on the proximity of the immune cells to the membrane. A more sophisticated approach might be to model the cells secreting proteases which diffuse freely and erode the membrane when they meet it. Whilst this is perhaps better fitted to the true processes that degrade the membrane, it is still not clear which cells are responsible for mediating this process. One potential route would thus be to consider the effects of having different immune cells responsible for this process within the modeling framework.

5.1.6. Waxing and Waning

One of the phenomena not accounted for in our description is the honeymoon period (Akirav et al., 2008). This is a phase lasting, in some cases, up to and beyond a year, typically immediately following diagnosis, in which the requirement for exogenous insulin is diminished. Though understanding of the mechanisms giving rise to this transient asymptomatic period is poor, one notable theoretical study demonstrated that this can occur naturally (in rare cases) when taking into account the role of endoplasmic reticulum stress and subsequent beta cell apoptosis (Jaberi-Douraki et al., 2014c). This study highlighted the possibility that elevations in the maximal unfolded protein response due to metabolic therapies could benefit patients by helping to meet metabolic demand, but despite this, they could not ultimately maintain a sufficient population of healthy beta cells to regulate blood glucose levels.

It is not clear whether the honeymoon period is reflective of events occurring at the level of individual islets, whether it is only present at the systemic level and whether it arises as a result of some innate periodicity present in either the immune system or in the beta cells. Modeling studies involving multiple islets may help to address some of these questions.

5.1.7. Therapeutic Intervention

Once a realistic geometry, motivated by real data, has been constructed, and the model sufficiently developed, an ultimate goal would be to test therapeutic interventions that target the pancreas specifically, such as the transplantation/implantation of encapsulated islets or populations of beta cells (Rickels et al., 2005; Robertson, 2010). We have discussed one hypothesis involving the promotion of low avidity T cell clones (Amrani et al., 2000; Khadra et al., 2009). One potential method for doing this is through the intravenous injection of ex vivo-selected and ex vivo-expanded autologous regulatory T cells, such as CD4+ cells (Bluestone et al., 2015), another is to incite the same expansion of regulatory T cells in vivo through the intravenous injection of nanoparticles coated with pMHC (Sugarman et al., 2013). Our modeling framework could easily be extended to include the effect of infusion of either the cells themselves or the nanoparticles to test the spatial aspects of these interventions.

5.1.8. Cell-Cell Coupling

It is evident from the multiple sources of data that communication between immune cells is an important aspect of insulitis, though the exact processes by which these cells exchange information is unknown. One potential avenue of investigation is to assess specific cell interactions based on known immunological mechanisms (Huppa and Davis, 2003; Friedl et al., 2005). These can be incorporated into an agent-based model to see what the overall effects of these specific types of interaction are, with a view to identifying which signaling pathways are most likely to result in the behavior observed in patient data (McLennan et al., 2012).

5.1.9. Comparison with Animal Models

Animal models, and in particular, the NOD mouse model have been extensively used to unravel potential causes and treatment options for human T1D. That true parallels between the disease in rodents and man may be fewer than is ideal, and the fact that treatments that work in mice are not as efficacious in humans (In't Veld, 2014) has led to questions about the suitability of animal models for human disease. From a modeling perspective, this opens up interesting questions. In particular, is it that the underlying mechanisms of insulitis truly are different between different species, or are some processes universal, albeit with potential different time courses and parameters? The former suggests that in order to find a cure for human T1D, only focus on the human condition is likely to yield fruitful results. However, the latter option suggests that, if relations between mice and men for those preserved pathways can be found, that the models can still provide vital insights. An in silico modeling framework seems like a natural place to address such a question as it allows quantitative probing of such relationships.

5.1.10. Experimental Design

One of the critical factors hindering progress in understanding T1D is the lack of human pancreatic samples from which to extract data. Moreover, these provide data only at one specific timepoint, so that very little is known about the time course of insulitis in humans.

One role for mathematical modeling is to identify avenues to explore to aid our understanding of disease mechanisms. The present study has suggested that the balance between repair and degradation of the peri-islet basement membrane surrounding the islets is a critical factor which determines the overall rate of beta cell destruction. Presently, little is known about the composition of this membrane in humans, though recent studies have identified its components in both healthy and infiltrated islets in the NOD mouse (Korpos et al., 2013). We thus feel that gaining an improved understanding of the structure of the human peri-islet basement membrane will significantly our improve understanding of insulitis.

Results from the agent-based model also suggest that the chemoattractant produced by the beta cells, and the lymphocytes' response to it is also important for governing disease progression. As such, in vitro experiments to assess and quantify the rate of immune cell migration to candidate chemokines would be useful to explore how immune cells initially reach the islet.

Finally, little is understood about the interaction between CD20+ and CD8+ cells during insulitis. The co-localization of these immune cell types in pancreatic samples suggest an interaction between the two (Huppa and Davis, 2003), and as such, we speculate that this could give rise to the activation of T cells. However, experiments have yet to confirm this for the immune cells involved in human T1D. Moreover, the specific dynamics of this process need to be better understood before the true role of the CD20+ cells in human insulitis can be revealed.

Overall, we believe that the proposed spatio-temporal framework has great potential as a tool to investigate insulitis in humans, and by focussing attention on the pancreas, it maximizes the use of the currently available human tissue data. We have identified a number of potential avenues in which the framework could be developed and a number of questions that could be addressed in doing so and we hope to continue in this line of research.

Author Contributions

All authors contributed to the initial design and subsequent development of the research project and model. KW developed code to run simulations and analyzed results from them. All authors drafted, edited and approved the submitted manuscript.

Funding

This work was generously supported by the Wellcome Trust Institutional Strategic Support Award (WT105618MA). KT gratefully acknowledges the financial support of the EPSRC via grant EP/N014391/1. We are also pleased to acknowledge financial support from the European Unions Seventh Framework Programme PEVNET [FP7/2007-2013] under grant agreement number 261441 to NM. The participants of the PEVNET consortium are described at http://www.uta.fi/med/pevnet/publications.html. Additional support was from a JDRF Career Development Award (5-CDA-2014-221-A-N) to SR and project grant 15/0005156 from Diabetes UK (to NM and SR).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys.2016.00633/full#supplementary-material

Video S1. Typical example of agent-based model simulation using parameters in Table 1, with immune cell lifespans of 28 days and B = 5.

Video S2. Typical example of agent-based model simulation using parameters in Table 1, with immune cell lifespans of 28 days and B = 30.

Video S3. Typical example of agent-based model simulation with no peri-islet membrane repair, as indicated in Section 4.2 with immune cell lifespans of 28 days, and B = 5.

Video S4. Typical example of agent-based model simulation with no peri-islet membrane repair, as indicated in Section 4.2 with immune cell lifespans of 28 days, and B = 30.

Video S5. Typical example of agent-based model simulation with high rates of peri-islet membrane degradation, as indicated in Section 4.3 with immune cell lifespans of 28 days, and B = 5.

Video S6. Typical example of agent-based model simulation with high rates of peri-islet membrane degradation, as indicated in Section 4.3 with immune cell lifespans of 28 days, and B = 30.

Video S7. Typical example of agent-based model simulation using in the high chemokine regime, as indicated in Section 4.5 with immune cell lifespans of 28 days, and B = 5.

Video S8. Typical example of agent-based model simulation using in the high chemokine regime, as indicated in Section 4.5 with immune cell lifespans of 28 days, and B = 30.

References

(2015). Type 1 Diabetes in Adults: Diagnosis Management. Technical Report August, National Institute for Health Care Excellence.

Akirav, E., Kushner, J. A., and Herold, K. C. (2008). Beta-cell mass and type 1 diabetes: going, going, gone? Diabetes 57, 2883–2888. doi: 10.2337/db07-1817

PubMed Abstract | CrossRef Full Text | Google Scholar

American Diabetes Association (2006). Diagnosis and classification of diabetes mellitus. Diabetes Care 29, S43–S48. doi: 10.2337/dc10-S062

PubMed Abstract | CrossRef Full Text

Amrani, A., Verdaguer, J., Serra, P., Tafuro, S., Tan, R., and Santamaria, P. (2000). Progression of autoimmune diabetes driven by avidity maturation of a T-cell population. Nature 406, 739–742. doi: 10.1038/35021081

PubMed Abstract | CrossRef Full Text | Google Scholar

An, G., Mi, Q., Dutta-Moscato, J., and Vodovotz, Y. (2009). Agent-based models in translational systems biology. Wiley Interdiscip. Rev. Syst. Biol. Med. 1, 159–171. doi: 10.1002/wsbm.45

PubMed Abstract | CrossRef Full Text | Google Scholar

Atkinson, M. A. (2012). The pathogenesis and natural history of type 1 diabetes. Cold Spring Harbor Perspect. Med. 2, 1–18. doi: 10.1101/cshperspect.a007641

PubMed Abstract | CrossRef Full Text | Google Scholar

Babon, J. A. B., DeNicola, M. E., Blodgett, D. M., Crèvecoeur, I., Buttrick, T. S., Maehr, R., et al. (2016). Analysis of self-antigen specificity of islet-infiltrating t cells from human donors with type 1 diabetes. Nat. Med. 22, 1482–1487. doi: 10.1038/nm.4203

PubMed Abstract | CrossRef Full Text | Google Scholar

Barnett, A. H., Eff, C., Leslie, R. D., and Pyke, D. A. (1981). Diabetes in identical twins. A study of 200 pairs. Diabetologia 20, 87–93. doi: 10.1007/bf00262007

PubMed Abstract | CrossRef Full Text | Google Scholar

Bertram, R., and Sherman, A. (2004). A calcium-based phantom bursting model for pancreatic islets. Bullet. Math. Biol. 66, 1313–1344. doi: 10.1016/j.bulm.2003.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Bluestone, J. A., Buckner, J. H., Fitch, M., Gitelman, S. E., Gupta, S., Hellerstein, M. K., et al. (2015). Type 1 diabetes immunotherapy using polyclonal regulatory T cells. Sci. Trans. Med. 7:315ra189. doi: 10.1126/scitranslmed.aad4134

PubMed Abstract | CrossRef Full Text | Google Scholar

Bogdani, M. (2016). Thinking outside the cell: a key role for hyaluronan in the pathogenesis of human type 1 diabetes. Diabetes 65, 2015–2114. doi: 10.2337/db15-1750

PubMed Abstract | CrossRef Full Text | Google Scholar

Boitard, C. (2012). Pancreatic islet autoimmunity. Presse Med. 41, e636–e650. doi: 10.1016/j.lpm.2012.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Bottazzo, G. F., Dean, B. M., McNally, J. M., MacKay, E. H., Swift, P. G., and Gamble, D. R. (1985). In situ characterization of autoimmune phenomena and expression of HLA molecules in the pancreas in diabetic insulitis. New Engl. J. Med. 313, 353–360. doi: 10.1056/NEJM198508083130604

PubMed Abstract | CrossRef Full Text | Google Scholar

Brehm, M. A., Powers, A. C., Shultz, L. D., and Greiner, D. L. (2012). Advancing animal models of human type 1 diabetes by engraftment of functional human tissues in immunodeficient mice. Cold Spring Harb. Perspect. Med. 2:a007757. doi: 10.1101/cshperspect.a007757

PubMed Abstract | CrossRef Full Text | Google Scholar

Brereton, M. F., Iberl, M., Shimomura, K., Zhang, Q., Adriaenssens, A. E., Proks, P., et al. (2014). Reversible changes in pancreatic islet structure and function produced by elevated blood glucose. Nat. Commun. 5:4639. doi: 10.1038/ncomms5639

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruna, M., and Chapman, S. J. (2012). Excluded-volume effects in the diffusion of hard spheres. Phys. Rev. E Stat. Nonline. Soft Matter Phys. 85, 24–29. doi: 10.1103/physreve.85.011103

PubMed Abstract | CrossRef Full Text | Google Scholar

Campbell-Thompson, M., Wasserfall, C., Kaddis, J., Albanese-O'Neill, A., Staeva, T., Nierras, C., et al. (2012). Network for Pancreatic Organ Donors with Diabetes (nPOD): developing a tissue biobank for type 1 diabetes. Diabetes Metab. Res. Rev. 28, 608–617. doi: 10.1002/dmrr.2316

PubMed Abstract | CrossRef Full Text | Google Scholar

Castiglioni, P., Gerloni, M., Cortez-Gonzales, X., and Zanetti, M. (2005). CD8 T cell priming by B lymphocytes is CD4 help dependent. Eur. J. Immunol. 35, 1360–1370. doi: 10.1002/eji.200425530

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, N., Glazier, J. A., Izaguirre, J. A., and Alber, M. S. (2007). A parallel implementation of the Cellular Potts Model for simulation of cell-based morphogenesis. Comput. Phys. Commun. 176, 670–681. doi: 10.1016/j.cpc.2007.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Christen, U., and Von Herrath, M. G. (2004). IP-10 and type 1 diabetes: a question of time and location. Autoimmunity 37, 273–282. doi: 10.1080/08916930410001713124

PubMed Abstract | CrossRef Full Text | Google Scholar

Cnop, M., Hughes, S. J., Igoillo-Esteve, M., Hoppa, M. B., Sayyed, F., van de Laar, L., et al. (2010). The long lifespan and low turnover of human islet beta cells estimated by mathematical modelling of lipofuscin accumulation. Diabetologia 53, 321–330. doi: 10.1007/s00125-009-1562-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Cnop, M., Welsh, N., Jonas, J.-C., Jo, A., and Lenzen, S. (2005). Mechanisms of pancreatic beta-cell death in type 1 and type 2 diabetes: many differences, few similarities. Diabetes 54, S97–S107. doi: 10.2337/diabetes.54.suppl_2.S97

PubMed Abstract | CrossRef Full Text | Google Scholar

Coppieters, K., Dotta, F., Amirian, N., Campbell, P. D., Kay, T. W. H., Atkinson, M. A., et al. (2012). Demonstration of islet-autoreactive CD8 t cells in insulitic lesions from recent onset and long-term type 1 diabetes patients. J. Exp. Med. 209, 51–60. doi: 10.1084/jem.20111187

PubMed Abstract | CrossRef Full Text | Google Scholar

Daub, J. T., and Merks, R. M. H. (2013). A cell-based model of extracellular-matrix-guided endothelial cell migration during angiogenesis. Bullet. Math. Biol. 75, 1377–1399. doi: 10.1007/s11538-013-9826-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Dirice, E., Kahraman, S., Jiang, W., El Ouaamari, A., De Jesus, D. F., Teo, A. K. K., et al. (2014). Soluble factors secreted by T cells promote β-cell proliferation. Diabetes 63, 188–202. doi: 10.2337/db13-0204

PubMed Abstract | CrossRef Full Text | Google Scholar

Dotta, F., Censini, S., van Halteren, A. G. S., Marselli, L., Masini, M., Dionisi, S., et al. (2007). Coxsackie B4 virus infection of beta cells and natural killer cell insulitis in recent-onset type 1 diabetic patients. Proc. Natl. Acad. Sci. U.S.A. 104, 5115–5120. doi: 10.1073/pnas.0700442104

PubMed Abstract | CrossRef Full Text | Google Scholar

Eisenbarth, G. S. (1986). Type I diabetes mellitus. A chronic autoimmune disease. New Engl. J. Med. 314, 1360–1368. doi: 10.1056/NEJM198605223142106

PubMed Abstract | CrossRef Full Text | Google Scholar

Elliott, C. M., Stinner, B., and Venkataraman, C. (2012). Modelling cell motility and chemotaxis with evolving surface finite elements. J. R. Soc. Interface 9, 3027–3044. doi: 10.1098/rsif.2012.0276

PubMed Abstract | CrossRef Full Text | Google Scholar

Foulis, A., and Stewart, J. (1984). The pancreas in recent-onset Type 1 (insulin-dependent) diabetes mellitus: insulin content of islets, insulitis and associated changes in the exocrine acinar tissue. Diabetologia 26, 456–461. doi: 10.1007/BF00262221

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedl, P., den Boer, A. T., and Gunzer, M. (2005). Tuning immune responses: diversity and adaptation of the immunological synapse. Nat. Rev. Immunol. 5, 532–545. doi: 10.1038/nri1647

PubMed Abstract | CrossRef Full Text | Google Scholar

Gepts, W. (1965). Pathologic anatomy of the pancreas in juvenile diabetes mellitus. Diabetes 14, 619–633. doi: 10.2337/diab.14.10.619

PubMed Abstract | CrossRef Full Text | Google Scholar

Graham, K. L., Krishnamurthy, B., Fynch, S., Ayala-Perez, R., Slattery, R. M., Santamaria, P., et al. (2012). Intra-islet proliferation of cytotoxic T lymphocytes contributes to insulitis progression. Eur. J. Immunol. 42, 1717–1722. doi: 10.1002/eji.201242435

PubMed Abstract | CrossRef Full Text | Google Scholar

Green, D. R., Droin, N., and Pinkoski, M. (2003). Activation-induced cell death in T cells. Immunol. Rev. 193, 70–81. doi: 10.1034/j.1600-065X.2003.00051.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, B., Serra, P., Yamanouchi, J., Amrani, A., Elliott, J. F., Dickie, P., et al. (2016). Developmental control of CD8+ t cell–avidity maturation in autoimmune diabetes. J. Clin. Invest. 115, 1879–1887. doi: 10.1172/JCI24219

PubMed Abstract | CrossRef Full Text | Google Scholar

Harvey, D. G., Fletcher, A. G., Osborne, J. M., and Pitt-Francis, J. (2015). A parallel implementation of an off-lattice individual-based model of multicellular populations. Comput. Phys. Commun. 192, 130–137. doi: 10.1016/j.cpc.2015.03.005

CrossRef Full Text | Google Scholar

Holcombe, M., Adra, S., Bicak, M., Chin, S., Coakley, S., Graham, A. I., et al. (2012). Modelling complex biological systems using an agent-based approach. Integr. Biol. 4, 53–64. doi: 10.1039/C1IB00042J

PubMed Abstract | CrossRef Full Text | Google Scholar

Huppa, J. B., and Davis, M. M. (2003). T-cell-antigen recognition and the immunological synapse. Nat. Rev. Immunol. 3, 973–983. doi: 10.1038/nri1245

PubMed Abstract | CrossRef Full Text | Google Scholar

In't Veld, P. (2014). Insulitis in human type 1 diabetes: a comparison between patients and animal models. Semin. Immunopathol. 36, 569–579. doi: 10.1007/s00281-014-0438-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Inzucchi, S. E. (2012). Diagnosis of diabetes. New Engl. J. Med. 367, 542–550. doi: 10.1056/NEJMcp1103643

PubMed Abstract | CrossRef Full Text | Google Scholar

Itoh, N., Hanafusa, T., Miyazaki, A., Miyagawa, J., Yamagata, K., Yamamoto, K., et al. (1993). Mononuclear cell infiltration and its relation to the expression of major histocompatibility complex antigens and adhesion molecules in pancreas biopsy specimens from newly diagnosed insulin-dependent diabetes mellitus patients. J. Clin. Invest. 92, 2313–2322. doi: 10.1172/JCI116835

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaberi-Douraki, M., Liu, S. W., Pietropaolo, M., and Khadra, A. (2014a). Autoimmune responses in T1DM: quantitative methods to understand onset, progression and prevention of disease. Pediatr. Diabetes 15, 162–174. doi: 10.1111/pedi.12148

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaberi-Douraki, M., Pietropaolo, M., and Khadra, A. (2014b). Predictive models of type 1 diabetes progression: understanding T-cell cycles and their implications on autoantibody release. PLoS ONE 9:e93326. doi: 10.1371/journal.pone.0093326

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaberi-Douraki, M., Pietropaolo, M., and Khadra, A. (2015). Continuum model of T-cell avidity: understanding autoreactive and regulatory T-cell responses in type 1 diabetes. J. Theor. Biol. 383, 93–105. doi: 10.1016/j.jtbi.2015.07.032

PubMed Abstract | CrossRef Full Text

Jaberi-Douraki, M., Schnell, S., Pietropaolo, M., and Khadra, A. (2014c). Unraveling the contribution of pancreatic beta-cell suicide in autoimmune type 1 diabetes. J. Theor. Biol. 375, 77–87. doi: 10.1016/j.jtbi.2014.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Janeway, C. A. Jr. Bottomly, K. I. M., Horowitz, J. A. Y., Kaye, J., Jones, B., and Tite, J. (1985). Modes of cell: cell communication in the immune system. J. Immunol. 135, 739s–742s.

PubMed Abstract | Google Scholar

Jones, A. G., and Hattersley, A. T. (2013). The clinical utility of C-peptide measurement in the care of patients with diabetes. Diabetic Med. 30, 803–817. doi: 10.1111/dme.12159

PubMed Abstract | CrossRef Full Text | Google Scholar

Kachapati, K., Adams, D., Bednar, K., and Ridgway, W. M. (2012). The non-obese diabetic (NOD) mouse as a model of human type 1 diabetes. Methods Mol. Biol. 933, 3–16. doi: 10.1007/978-1-62703-068-7_1

PubMed Abstract | CrossRef Full Text | Google Scholar

Khadra, A., Pietropaolo, M., Nepom, G. T., and Sherman, A. (2011). Investigating the role of T-cell avidity and killing efficacy in relation to type 1 diabetes prediction. PLoS ONE 6:e14796. doi: 10.1371/journal.pone.0014796

PubMed Abstract | CrossRef Full Text | Google Scholar

Khadra, A., Santamaria, P., and Edelstein-Keshet, L. (2009). The role of low avidity T cells in the protection against type 1 diabetes: a modeling investigation. J. Theor. Biol. 256, 126–141. doi: 10.1016/j.jtbi.2008.09.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, A., Miller, K., Jo, J., Kilimnik, G., and Wojcik, P. (2009). Islet architecture. Islets 1, 129–136. doi: 10.4161/isl.1.2.9480

PubMed Abstract | CrossRef Full Text | Google Scholar

Klöppel, G., Löhr, M., Habich, K., Oberholzer, M., and Heitz, P. U. (1985). Islet pathology and the pathogenesis of type 1 and type 2 diabetes mellitus revisited. Surv. Synth. Pathol. Res. 4, 110–125.

PubMed Abstract | Google Scholar

Knip, M., Veijola, R., Virtanen, S. M., Hyo, H., Vaarala, O., and Åkerblom, H. K. (2005). Environmental triggers and determinants of type 1 diabetes. Diabetes 54, S125–S136. doi: 10.2337/diabetes.54.suppl_2.S125

PubMed Abstract | CrossRef Full Text | Google Scholar

Korpos, É., Kadri, N., Kappelhoff, R., Wegner, J., Overall, C. M., Weber, E., et al. (2013). The peri-islet basement membrane, a barrier to infiltrating leukocytes in type 1 diabetes in mouse and human. Diabetes 62, 531–542. doi: 10.2337/db12-0432

PubMed Abstract | CrossRef Full Text | Google Scholar

La Torre, D., and Lernmark, A. (2012). Immunology of beta cell destruction. Adv. Exp. Med. Biol. 771, 537–583. doi: 10.1007/978-90-481-3271-3_24

CrossRef Full Text | Google Scholar

Lally, F. J., and Bone, A. J. (2003). “Animal models of type 1 diabetes,” in Textbook of Diabetes, 3rd Edn., eds J. C. Pickup and G. Williams (Oxford: Blackwell Science), 19.1–19.17.

Google Scholar

Lecompte, P. (1958). Insulitis in early juvenile diabetes. AMA Arch. Pathol. 66, 450–457.

PubMed Abstract | Google Scholar

Leete, P., Willcox, A., Krogvold, L., Dahl-Jørgensen, K., Foulis, A. K., Richardson, S. J., et al. (2016). Differential insulitic profiles determine the extent of beta cell destruction and the age at onset of type 1 diabetes. Diabetes 65, 1362–1369. doi: 10.2337/db15-1615

PubMed Abstract | CrossRef Full Text

Levine, H., Rappel, W.-J., and Cohen, I. (2000). Self-organization in systems of self-propelled particles. Phys. Rev. E 63:017101. doi: 10.1103/PhysRevE.63.017101

PubMed Abstract | CrossRef Full Text | Google Scholar

Lo, S. S., Tun, R. Y., Hawa, M., and Leslie, R. D. (1991). Studies of diabetic twins. Diabetes Metab. Res. Rev. 7, 223–238. doi: 10.1002/dmr.5610070403

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahaffy, J. M., and Edelstein-Keshet, L. (2007). Modeling cyclic waves of circulating t cells in autoimmune diabetes. SIAM J. Appl. Math. 67, 915–937. doi: 10.1137/060661144

CrossRef Full Text | Google Scholar

Maini, P. K., and Baker, R. E. (2011). Modelling collective cell motion in biology. Adv. Appl. Math. 47, 688–709. doi: 10.1007/978-3-319-06923-4_1

CrossRef Full Text | Google Scholar

Marée, A., Grieneisen, V., and Hogeweg, P. (2007). “Single-cell based models in biology and medicine,” in The Cellular Potts Model and Biophysical Properties of Cells, Tissues and Morphogenesis (Basel: Birkhäuser Verlag), 107–136.

PubMed Abstract

Marée, A. F. M., Kublik, R., Finegood, D. T., and Edelstein-Keshet, L. (2006). Modelling the onset of Type 1 diabetes: can impaired macrophage phagocytosis make the difference between health and disease? Philos. Trans. Ro. Soc. A 364 1267–1282. doi: 10.1098/rsta.2006.1769

PubMed Abstract | CrossRef Full Text

McLennan, R., Dyson, L., Prather, K. W., Morrison, J. A., Baker, R. E., Maini, P. K., et al. (2012). Multiscale mechanisms of cell migration during development: theory and experiment. Development 139, 2935–2944. doi: 10.1242/dev.081471

PubMed Abstract | CrossRef Full Text

Merks, R. M. H., Brodsky, S. V., Goligorsky, M. S., Newman, S. A., and Glazier, J. A. (2006). Cell elongation is key to in silico replication of in vitro vasculogenesis and subsequent remodeling. Dev. Biol. 289, 44–54. doi: 10.1016/j.ydbio.2005.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Middleton, A. M., Fleck, C., and Grima, R. (2014). A continuum approximation to an off-lattice individual-cell based model of cell migration and adhesion. J. Theor. Biol. 359, 220–232. doi: 10.1016/j.jtbi.2014.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirams, G., Arthurs, C., Bernabeu, M., Bordas, R., Cooper, J., Corrias, A., et al. (2013). Chaste: an open source C++ library for computational physiology and biology. PLOS Comput. Biol. 9:e1002970. doi: 10.1371/journal.pcbi.1002970

PubMed Abstract | CrossRef Full Text | Google Scholar

Morgan, N. G., Leete, P., Foulis, A. K., and Richardson, S. J. (2014). Islet inflammation in human type 1 diabetes mellitus. IUBMB Life 66, 723–734. doi: 10.1002/iub.1330

PubMed Abstract | CrossRef Full Text | Google Scholar

Palsson, E., and Othmer, H. G. (2000). A model for individual and collective cell movement in Dictyostelium discoideum. Proc. Natl. Acad. Sci. U.S.A. 97, 10448–10453. doi: 10.1073/pnas.97.19.10448

PubMed Abstract | CrossRef Full Text | Google Scholar

Plank, M. J., and Simpson, M. J. (2012). Models of collective cell behaviour with crowding effects: comparing lattice-based and lattice-free approaches. J. R. Soc. Interface 9, 2983–2996. doi: 10.1098/rsif.2012.0319

PubMed Abstract | CrossRef Full Text | Google Scholar

Pugliese, A. (2014). Advances in the etiology and mechanisms of type 1 diabetes. Discov. Med. 18, 141–150.

PubMed Abstract | Google Scholar

Pugliese, A., Vendrame, F., Reijonen, H., Atkinson, M. A., Campbell-Thompson, M., and Burke, G. W. (2014). New insight on human type 1 diabetes Biology: nPOD and nPOD-transplantation. Curr. Diabetes Rep. 14:530. doi: 10.1007/s11892-014-0530-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Redondo, M. J., Yu, L., Hawa, M., Mackenzie, T., Pyke, D. A., Eisenbarth, G. S., et al. (2001). Heterogeneity of type I diabetes: analysis of monozygotic twins in Great Britain and the United States. Diabetologia 44, 354–62. doi: 10.1007/s001250051626

PubMed Abstract | CrossRef Full Text | Google Scholar

Reed, J. C., and Herold, K. C. (2015). Thinking bedside at the bench: the NOD mouse model of T1DM. Nat. Rev. Endocrinol. 11, 308–314. doi: 10.1038/nrendo.2014.236

PubMed Abstract | CrossRef Full Text | Google Scholar

Reijonen, H., Mallone, R., Heninger, A.-K., Laughlin, E. M., Kochik, S. A., Falk, B., et al. (2004). GAD65-specific CD4 t-cells with high antigen avidity are prevalent in peripheral blood of patients with type 1 diabetes. Diabetes 53, 1987–1994. doi: 10.2337/diabetes.53.8.1987

PubMed Abstract | CrossRef Full Text | Google Scholar

Richardson, S. J., Morgan, N. G., and Foulis, A. K. (2014). Pancreatic pathology in type 1 diabetes mellitus. Endocr. Pathol. 25, 80–92. doi: 10.1007/s12022-014-9297-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Richardson, S. J., Willcox, A., Bone, A. J., Morgan, N. G., and Foulis, A. K. (2011). Immunopathology of the human pancreas in type-I diabetes. Semin. Immunopathol. 33, 9–21. doi: 10.1007/s00281-010-0205-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Rickels, M. R., Schutta, M. H., Markmann, J. F., Barker, C. F., Naji, A., and Teff, K. L. (2005). β-Cell function following human islet transplantation for type 1 diabetes. Diabetes 54, 100–106. doi: 10.2337/diabetes.54.1.100

PubMed Abstract | CrossRef Full Text | Google Scholar

Robertson, R. P. (2010). Islet transplantation a decade later and strategies for filling a half-full glass. Diabetes 59, 1285–1291. doi: 10.2337/db09-1846

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodríguez-Pinto, D. (2005). B cells as antigen presenting cells. Cell. Immunol. 238, 67–75. doi: 10.1016/j.cellimm.2006.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Roep, B. O., and Tree, T. I. M. (2014). Immune modulation in humans: implications for type 1 diabetes mellitus. Nat. Rev. Endocrinol. 10, 229–242. doi: 10.1038/nrendo.2014.2

PubMed Abstract | CrossRef Full Text | Google Scholar

Saisho, Y., Butler, A. E., Manesso, E., Elashoff, D., Rizza, R. A., and Butler, P. C. (2013). Beta-cell mass and turnover in humans. Diabetes Care 36, 111–117. doi: 10.2337/dc12-0421

PubMed Abstract | CrossRef Full Text | Google Scholar

Schanze, T. (1995). Sinc interpolation of discrete periodic signals. IEEE Trans. Signal Process. 43, 1502–1503. doi: 10.1109/78.388863

CrossRef Full Text

Scianna, M. (2015). An extended cellular potts model analyzing a wound healing assay. Comput. Biol. Med. 62, 33–54. doi: 10.1016/j.compbiomed.2015.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Somoza, N., Vargas, F., Roura-Mir, C., Vives-Pi, M., Fernández-Figueras, M. T., Ariza, A., et al. (1994). Pancreas in recent onset insulin-dependent diabetes mellitus. Changes in HLA, adhesion molecules and autoantigens, restricted T cell receptor V beta usage, and cytokine profile. J. Immunol. 153, 1360–1377.

PubMed Abstract | Google Scholar

Song, C., Wang, P., and Makse, H. A. (2008). A phase diagram for jammed matter. Nature 453, 629–632. doi: 10.1038/nature06981

PubMed Abstract | CrossRef Full Text | Google Scholar

Sprent, J. (1993). Lifespans of naive, memory and effector lymphocytes. Curr. Opin. Immunol. 5, 433–438. doi: 10.1016/0952-7915(93)90065-Z

PubMed Abstract | CrossRef Full Text | Google Scholar

Standifer, N., Burwell, E. A., Gersuk, V. H., Greenbaum, C. J., and Nepom, G. T. (2009). Changes in autoreactive t cell avidity during type 1 diabetes development. Clin. Immunol. 132, 312–320. doi: 10.1016/j.clim.2009.04.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Stein, J. V., and Nombela-Arrieta, C. (2005). Chemokine control of lymphocyte trafficking: a general overview. Immunology 116, 1–12. doi: 10.1111/j.1365-2567.2005.02183.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugarman, J., Tsai, S., Santamaria, P., and Khadra, A. (2013). Quantifying the importance of pmhc valency, total pmhc dose and frequency on nanoparticle therapeutic efficacy. Immunol. Cell Biol. 91, 350–359. doi: 10.1038/icb.2013.9

PubMed Abstract | CrossRef Full Text | Google Scholar

Swat, M., Thomas, G. L., Belmonte, J. M., Shirinifard, A., Hmeljak, D., and Glazier, J. A. (2012). Multi-scale modeling of tissues using CompuCell3D. Comput. Methods Cell Biol. Methods Cell Biol. 110, 325–366. doi: 10.1016/B978-0-12-388403-9.00013-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Tapia, J. J., and D'souza, R. M. (2011). Parallelizing the cellular potts model on graphics processing units. Comput. Phys. Commun. 182, 857–865. doi: 10.1016/j.cpc.2010.12.011

CrossRef Full Text | Google Scholar

Taplin, C. E., and Barker, J. M. (2008). Autoantibodies in type 1 diabetes. Autoimmunity 41, 11–18. doi: 10.1080/08916930701619169

PubMed Abstract | CrossRef Full Text | Google Scholar

Tijskens, E., Ramon, H., and De Baerdemaeker, J. (2003). Discrete element modelling for process simulation in agriculture. J. Sound Vibr. 266, 493–514. doi: 10.1016/S0022-460X(03)00581-9

CrossRef Full Text | Google Scholar

Trudeau, J. D., Kelly-smith, C., Verchere, C. B., Elliott, J. F., Dutz, J. P., Finegood, D. T., et al. (2003). Prediction of spontaneous autoimmune diabetes in nod mice by quantification of autoreactive t cells in peripheral blood. J. Clin. Invest. 111, 217–223. doi: 10.1172/JCI200316409

PubMed Abstract | CrossRef Full Text | Google Scholar

Walker, J. N., Johnson, P. R. V., Shigeto, M., Hughes, S. J., Clark, A., and Rorsman, P. (2011). Glucose-responsive beta cells in islets isolated from a patient with long-standing type 1 diabetes mellitus. Diabetologia 54, 200–202. doi: 10.1007/s00125-010-1930-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Abbasi, F., Ornatsky, O., Cole, K. D., Misakian, M., Gaigalas, A. K., et al. (2012). Human CD4+ lymphocytes for antigen quantification: characterization using conventional flow cytometry and mass cytometry. Cytometry A 81A, 567–575. doi: 10.1002/cyto.a.22060

CrossRef Full Text | Google Scholar

Willcox, A., Richardson, S. J., Bone, A. J., Foulis, A. K., and Morgan, N. G. (2009). Analysis of islet inflammation in human type 1 diabetes. Clin. Exp. Immunol. 155, 173–181. doi: 10.1111/j.1365-2249.2008.03860.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Willcox, A., Richardson, S. J., Bone, A. J., Foulis, A. K., and Morgan, N. G. (2010). Evidence of increased islet cell proliferation in patients with recent-onset type 1 diabetes. Diabetologia 53, 2020–2028. doi: 10.1007/s00125-010-1817-6

PubMed Abstract | CrossRef Full Text | Google Scholar