## Abstract

A wave of apoptosis (programmed cell death) occurs normally in pancreatic β-cells of newborn mice. We previously showed that macrophages from non-obese diabetic (NOD) mice become activated more slowly and engulf apoptotic cells at a lower rate than macrophages from control (Balb/c) mice. It has been hypothesized that this low clearance could result in secondary necrosis, escalating inflammation and self-antigen presentation that later triggers autoimmune, Type 1 diabetes (T1D). We here investigate whether this hypothesis could offer a reasonable and parsimonious explanation for onset of T1D in NOD mice. We quantify variants of the *Copenhagen model* (Freiesleben De Blasio *et al*. 1999 *Diabetes* **48**, 1677), based on parameters from NOD and Balb/c experimental data. We show that the original Copenhagen model fails to explain observed phenomena within a reasonable range of parameter values, predicting an unrealistic all-or-none disease occurrence for both strains. However, if we take into account that, in general, activated macrophages produce harmful cytokines *only* when engulfing necrotic (but not apoptotic) cells, then the revised model becomes qualitatively and quantitatively reasonable. Further, we show that known differences between NOD and Balb/c mouse macrophage kinetics are large enough to account for the fact that an apoptotic wave can trigger escalating inflammatory response in NOD, but not Balb/c mice. In Balb/c mice, macrophages clear the apoptotic wave so efficiently, that chronic inflammation is prevented.

## 1. Introduction

In this paper, we explore the initial events underlying Type 1 diabetes (T1D), an autoimmune disease in which the destruction of insulin producing β-cells is due to an abnormal immune response (Beyan *et al*. 2003). The major destructive mediators are specific clones of cytotoxic T-cells that invade the pancreatic islets of Langerhans (Schranz & Lernmark 1998; Höglund *et al*. 1999; Devendra *et al*. 2004). However, many other types of immune cells and chemicals are implicated in setting the stage for T1D. Among these are macrophages, the professional phagocytes of the immune system and antigen-presenting cells (e.g. dendritic cells).

Overt diabetes, characterized by a dramatic increase in the blood sugar level, is preceded by insulitis, an inflammation of the pancreas. The event triggering inflammation could be viral infection, toxins, dietary factors (Schranz & Lernmark 1998; Beyan *et al*. 2003), or other causes. Numerous genetic and environmental factors correlate with T1D. In humans, this is difficult to study non-invasively (Schranz & Lernmark 1998; Thomas & Kay 2000). Animal models for T1D include the non-obese diabetic (NOD) mouse, used extensively to investigate the initiation and development of the disease. In NOD mice, immune cells infiltrate the pancreas and accumulate around the islets starting as early as two weeks after birth (Fujino-Kurihara *et al*. 1985). Macrophages and dendritic cells appear early (2–5 weeks), followed by CD4 and CD8 T-cells (6–10 weeks) (Dahlén *et al*. 1998). Full-blown diabetes ensues at about 12–15 weeks of age (Höglund *et al*. 1999) in female NOD mice, when 95% of the β-cells have been destroyed (Luppi & Trucco 1999).

In this paper, we concentrate on macrophage function at early stages of the disease. An important role of macrophages is to remove apoptotic cells and extracellular debris. This is of relevance to T1D, since death of β-cells occurs by apoptosis (Thomas & Kay 2000; Eizirik & Mandrup-Poulsen 2001), a form of cell death generally not associated with an inflammatory response. (In contrast, unregulated death, necrosis, is generally pro-inflammatory.) If apoptotic cells are not removed promptly, they become necrotic. Thus, defective phagocytosis by macrophages could lead to the build-up of necrotic material and escalating inflammation.

Our investigation is driven by several pre-existing experimental and theoretical studies. (i) It has been shown that a high rate of apoptosis of β-cells occurs as a natural part of early development. Trudeau *et al*. (2000) describe a wave of β-cell apoptosis in neonatal rodents, including NOD and Balb/c mice. In all cases, the wave peaks around two weeks of age. It has been estimated that over its course, up to 60% of the pre-existing β-cells die, with a peak rate of 9% per day. This study led to the hypothesis that the apoptotic wave could serve as a stimulus initiating inflammation (Trudeau *et al*. 2000; Mathis *et al*. 2001). (A possible loss of α-cells (Kaung 1994), which would increase the load of apoptotic cells to be cleared, will not be considered here; we will show that triggering an escalating inflammatory response does not depend on the actual size of the apoptotic wave, as long as it is above a certain threshold.) (ii) Further experiments showed that macrophages from NOD mice have defective phagocytosis, compared to healthy controls (Balb/c mice) (O'Brien *et al*. 2002*a*,*b*; Beyan *et al*. 2003), suggesting that inadequate macrophage phagocytosis might cause amplification of the initial inflammation. (iii) In a joint experimental–modelling study, we recently quantified details of phagocytic engulfment and digestion *in vitro* by macrophages from the two mouse strains (Marée *et al*. 2005). Data analysis revealed that Balb/c macrophages have a faster basal engulfment rate, and further undergo ‘activation’—accelerate this rate—upon encountering apoptotic cells. In NOD mouse macrophages, basal engulfment rates are low and do not change. The latter work has provided detailed quantitative comparisons of macrophage function that we now apply to parametrization of models in this paper, and comparison of strains.

In a pre-existing work, Freiesleben De Blasio *et al*. (1999) modelled interactions of macrophages, antigen and T-cells. They showed that their system, henceforth known as the *Copenhagen model*, can become unstable, resulting in autoimmunity. The main aim of the Copenhagen model was to postulate a qualitative hypothesis, rather than to make a quantitative prediction. Here we extend this idea to quantitative analysis based on mouse parameters. Our goal is not to model T1D in its full complexity, but rather to give a parsimonious explanation of its initiation, based on a limited set of well-established quantitative differences between NOD and Balb/c mice. We use rate constants determined by Marée *et al*. (2005) for macrophage engulfment kinetics to answer the following questions.

Can the differences in macrophage phagocytosis function in NOD versus Balb/c mice (alone, or in combination with other factors) account for the distinct fates of these two strains, i.e. possible initiation of autoimmunity in NOD but not in Balb/c mice?

Can the wave of β-cell death associated with normal development in all mice be a triggering stimulus that initiates the inflammation in NOD mice?

Based on the fact that not all NOD mice develop diabetes (or even chronic inflammation), any model describing early dynamics should minimally have a stable healthy rest state for *both* NOD and Balb/c mice. Further, the possibility of chronic inflammation in NOD mice can be interpreted in the model by a saddle point and a separatrix that separates regimes attracted to the healthy state from those that lead to chronic inflammation. A triggering stimulus would then be an impulse that carries the system into the basin of attraction of the inflamed state. Thus, our questions above can be reworded as follows: is such a dynamic configuration possible (for NOD but not Balb/c mice) within reasonable, biologically relevant ranges of the rate constants? We show that the answers to our questions are positive, but that the inflammatory threshold as proposed in the original Copenhagen model is an incomplete description, that cannot lead to this behaviour *within biologically relevant parameter ranges*. Only by incorporating other realistic features, such as necrosis, are both the quantitative and the qualitative aspect of this model accurate.

Here we are not concerned with dynamics at later stages, in which the adaptive immune system is implicated. Further, we make no claims about the corresponding human systems. However, it is to be hoped that once the development of diabetes in NOD mice is quantitatively well understood, we would be better situated to understand the dynamics in humans as well.

## 2. A basic model for macrophage-induced inflammation

We consider first a simple modification of the Copenhagen model, defined in figure 1, table 1 and equations (2.1*a*)–(2.1*c*). Denoting by *M* and *M*_{a} the density of resting and activated macrophages in the islets, and *B*_{a} the density of apoptotic cells, we formulate the following mass-action equations of balance (see figure 1 and table 1 for all definitions):(2.1a)

(2.1b)

(2.1c)

Essential features of the model that carry over from the original version by Freiesleben De Blasio *et al*. (1999) are the influx of macrophages (enhanced by recruitment); their interconversion between resting and ‘activated’ states; and their dual role in clearing apoptotic cells while simultaneously producing further damage. Our modifications of that model are as follows: (i) explicit consideration of β-cells in place of generalized ‘antigen’. (We can thereby incorporate the normal developmental apoptotic wave of *β* cells, *W*(*t*), quantified by Trudeau *et al*. (2000).) (ii) Inclusion of apoptotic cell removal by both resting and activated macrophages, as quantified by Marée *et al*. (2005). (iii) Inclusion of terms with parameters *e*_{i} to represent crowding effects (i.e. reduced entry and/or increased efflux of macrophages from the tissue at high densities.) These terms address unlimited growth that occurs in the original Copenhagen model. (We generally expect that these terms are significant only at high densities.) Further, we here interpret the activated state of macrophages to mean a state of accelerated phagocytosis (see Marée *et al*. 2005). In Balb/c mice, a basal rate of macrophage engulfment, *f*_{1}, is (as a consequence of engulfment) accelerated to a higher level, *f*_{2}. We consider this acceleration to be the ‘activation’ step (in contrast with other meanings in the literature, associated with differentiation of macrophages into professional antigen-presenting cells). We also assume that this activated state leads to secretion of cytokines that may harm β-cells and induce increased apoptosis. Our model reduces to the original Copenhagen model for *e*_{1}=*e*_{2}=0, *W*(*t*)=0 and *f*_{1}=0.

We analyse variants of this model in the following sections and show that, while, in principle, it can lead to dynamics described in our introduction, it *de facto* cannot do so in relevant parameter regimes. This analysis leads to a modification where specific influence of necrotic cells is included in the inflammatory cascade.

## 3. Simplest reduced model

We first comment on a stripped-down version of the model from which insight is obtained. In the case that macrophages neither enter nor leave the tissue (*J*=0, *c*=0, *b*=0), and are not affected by overcrowding (*e*_{1}=*e*_{2}=0), the model becomes easy to understand. Further omitting the stimulus (*W*(*t*)=0) and the role of resting macrophages in phagocytosis (*f*_{1}=0) leaves the simplest ‘interesting’ system:(3.1a)(3.1b)(3.1c)Then the total density of macrophages, , is constant, and consequently the model can be reduced to two equations ( and as above), and analysed fully by elementary phase-plane methods. A healthy rest state occurs at (*M*_{a}, *B*_{a})=(0, 0), but is unstable when a state of chronic inflammation at exists (i.e. when ). The condition for existence and stability of the inflamed state is(3.2)This is analogous to a threshold condition in susceptible–infectious–susceptible (SIS) epidemiology models. *r*_{0} is a ‘basic reproductive parameter for inflammation’, i.e. the average number of secondary activations caused by one activated macrophage (see table 2). Note that this parameter is equivalent to the ‘amount of secondary antigen produced by the primary β-cell destruction’, denoted *f*_{0} in Freiesleben De Blasio *et al*. (1999). In §4, we show a similar threshold condition in our basic model (2.1*a*)–(2.1*c*). At low macrophage density, only the healthy state can exist; at high density, relative to a tendency to transmit and respond to an ‘activation signal’ (represented by a parameter grouping), inflammation will occur. In SIS models, transmission of infection similarly depends on the susceptible pool size and parameters governing the transmission of infection.

## 4. Analysis of the basic model

We now consider the analysis of the model given by equations (2.1*a*)–(2.1*c*) in the case *W*(*t*)=0. Further, we assume (following Marée *et al*. (2005)) that *g*=*f*_{1}, which just means that resting macrophages become activated during the first engulfment of an apoptotic cell. In the case of negligible crowding effects of the resting macrophages, *e*_{1}=0, the system of equations (2.1*a*)–(2.1*c*) has a healthy rest state given by (*M*, *M*_{a}, *B*_{a})=(*J*/*c*, 0, 0), with a basal level of resting macrophages, *J*/*c*, and no apoptotic cells. When 0<*e*_{1}≪1, the healthy rest state, , can be approximated by terms in its series expansion (about *e*_{1}=0), , so that the approximation *M*≈*J*/*c* holds, provided *e*_{1}≪*c*^{2}/*J* (generally true, see table 1 and §5).

Analogous to *r*_{0} defined by equation (3.2), we can also define a threshold condition for the basic model that delineates when a low level of inflammation would be amplified. Close to (*J*/*c*, 0, 0), inflammation grows if , i.e. *gMB*_{a}−*kM*_{a}>0 or simply *gMB*_{a}/*kM*_{a}>1. Accordingly, define *r*_{0}=*gMB*_{a}/(*kM*_{a}) about the rest state. Then, using *M*≈*J*/*c*, , we simplify to obtain the condition for accelerating inflammation:(4.1)Freiesleben De Blasio *et al*. (1999) omit engulfment by resting macrophages (*f*_{1}=0), and scale their model, so that *J*=*d*=1. In that case, corresponds exactly to their equivalent threshold parameter.

The definition and interpretation of *r*_{0} is closely linked to stability of the healthy rest state. When *e*_{1} is small so that the approximation of the resting macrophage density by *J*/*c* holds, full stability analysis reveals that the healthy state is stable whenever(4.2)This is equivalent to *r*_{0}<1, i.e. complementary to the condition for growing inflammation (4.1).

It is cumbersome to find the non-trivial steady states of equations (2.1*a*)–(2.1*c*) because of the multiple nonlinearities. (Note that even when *e*_{1}≪*c*^{2}/*J*, the second term in the series for *M* can only be ignored close to the healthy rest state.) We therefore explore a related system of equations in which the crowding terms are absent, i.e. *e*_{1}=*e*_{2}=0. This system allows for at most one non-trivial positive equilibrium. Therefore, whenever the healthy equilibrium is stable, this equilibrium has to be a saddle point, which consequently separates the healthy regime from the inflamed regime. The simplification eliminates a finite high-density steady state (i.e. the inflammation becomes unbounded), but only mildly affects the existence and location of the saddle point. This ‘separating saddle point’ satisfies the equations: , , , and can be found explicitly, resulting in(4.3)Combining the constraint *M*, *M*_{a}, *B*_{a}>0 (for biological relevance of (4.3)), with stability of the healthy state (condition (4.2)), leads to(4.4a)(4.4b)This expresses that both strains have a healthy rest state, but only NOD mice can develop chronic inflammation. It is evident that these inequalities are restrictive, and here resides the difficulty in turning this conceptual model to a quantitatively feasible description of NOD mouse dynamics.

## 5. Parameter estimates

To obtain quantitative comparisons, we want to use realistic parameter estimates. In this section, we show how estimates for the parameters given in table 1 were obtained based on data from the literature, and our own previous experiments. To avoid variations due to volume size, cell densities carry units of cells per ml.

### (a) Rates of phagocytosis of apoptotic cells by macrophages

Engulfment rates were obtained to good precision as output of *in vitro* experiments and data-fitting reported in Marée *et al*. (2005).

The rate of phagocytosis of apoptotic cells by resting macrophages (and simultaneous activation of those macrophages) is

*f*_{1}=*g*=2×10^{−5}(for Balb/c) and*f*_{1}=*g*=1×10^{−5}(for NOD) ml cell^{−1}d^{−1}.The rate of phagocytosis of apoptotic cells by activated macrophages is

*f*_{2}=5×10^{−5}(for Balb/c) and 1×10^{−5}(for NOD) ml cell^{−1}d^{−1}.Digestion of engulfed apoptotic cells takes place serially at rate

*k*_{d}≈25 d^{−1}.

### (b) Cell densities, fluxes and turnover rates

The following parameter values are given as ballpark estimates unless otherwise indicated.

A mouse pancreatic islet is about 150 μm diameter, and the volume is thus roughly 1.77×10

^{6}μm^{3}≈1.8×10^{−6}ml (Rosmalen*et al*. 2000). An islet contains approximately 500–1000 β-cells. The β-cell density is thus in the range of 4×10^{8}cells ml^{−1}. The apoptotic wave peaks around age 1.5 weeks, when the rate of cell loss is around 9% per day (Trudeau*et al*. 2000; O'Brien*et al*. 2002*a*). (No differences between NOD and Balb/c mice have been observed.) This rate was established indirectly, by subtracting β-cell replication, but not neogenesis (which could not be measured), from the growth of the β-cell mass. Assuming a rate of neogenesis around 1%, we estimate that at the peak of the apoptotic wave about (0.01+0.09)×4×10^{8}=4×10^{7}cells ml^{−1}enter apoptosis each day. According to O'Brien*et al*. (2002*b*), in female NOD mice during the wave, 0.57% of the β-cells are TUNEL positive, i.e. the apoptotic β-cell density peaks around 2.4×10^{6}cells ml^{−1}. In Balb/c mice, however, this density is only one-third as large, i.e. about 8×10^{5}cells ml^{−1}. Taken together, this implies that most apoptotic β-cells are removed within a couple of hours (but much less efficiently in NOD mice). We argue that at this peak, the density of apoptotic β-cells must be far higher than any later level, since at the early neonatal stage, there is only a low level of clearance of the apoptotic cells pending recruitment of macrophages. (Further, if the apoptosis rate stayed this high, the total β-cell mass would be depleted before the adaptive immune system responds, contradicting evidence for the central role of T-cells.)According to Van Furth

*et al*. (1973), the normal (uninflamed) influx of macrophages is*J*=48 000≈5×10^{4}cells ml^{−1}d^{−1}. These macrophages turn over in 4–15 days (Van Furth 1989; Paul 1993). Similar arguments in Van Furth & Diesselhoff-den Dulk (1984) suggest that the turnover rate is in the range*c*=0.07–0.25 d^{−1}. We approximate*c*≈0.1 d^{−1}. (This contrasts with a lower estimate of*c*=0.011 d^{−1}cited by Wigginton & Kirschner (2001).) Combining influx and turnover in the healthy state, . The macrophage density in inflamed tissue is of the order of magnitude of*M*_{inflamed}≈1×10^{7}cells ml^{−1}.Apoptotic cells undergo necrosis in 1–2 days (Van Nieuwenhuijze

*et al*. 2003), so the rate of necrotic degradation of apoptotic β-cells is*d*≈0.5–1 d^{−1}.In our previous experiments, we quantified successive steps of engulfment (and digestion) of apoptotic cells by activated macrophages. Interpreting activated macrophages as those with one or more engulfed apoptotic cells, we described the macrophage classes

*M*_{1}⇋*M*_{2}⇋*M*_{3}…, where the subscripts correspond to numbers of engulfed cells. Transitions to the right (engulfments of apoptotic cells, rate*f*_{2}*B*_{a}) and to the left (digestion, rate*k*_{d}) imply that, at quasi-steady state (QSS),*M*_{2}=*λM*_{1}and the number of macrophages in class*n*is*M*_{n}=*λM*_{n−1}, where*λ*=*f*_{2}*B*_{a}/*k*_{d}. (See details in Marée*et al*. (2005).) We find that . However, only those in class*M*_{1}can revert to a resting state. These make up a fraction*M*_{1}/*M*_{a}=(1−*λ*) of activated macrophages. An estimate for the deactivation rate is then . Consequently, when*f*_{2}*B*_{a}/*k*_{d}1, the deactivation rate*k*becomes much smaller than*k*_{d}≈25 d^{−1}; this is actually the case for both Balb/c and NOD mice during the whole apoptotic wave (i.e. when*B*_{a}>4×10^{5}for Balb/c, or*B*_{a}>1×10^{6}for NOD,*k*<5 d^{−1}). (This also means that the inactivation rate should depend on the number of apoptotic cells, with effectively no inactivation taking place when*B*_{a}densities are very high; here, however, we ignore this further complication.) Moreover, studies have shown that a complete digestion of engulfed apoptotic cells does not imply an instantaneous return to the resting state, so the above should be considered to be an upper bound. For example, Licht*et al*. (2001) and O'Brien*et al*. (2002*b*) found increases in the phagocytic abilities of macrophages, following uptake of apoptotic cells a day earlier. Evidence for the time-scale of macrophage inactivation after digestion of the apoptotic meal is controversial, but likely to be of the order of 24–48 h. Taken together, we estimate*k*≈0.4 d^{−1}.The parameter

*b*represents recruitment of resting macrophages into the tissue due to the presence of activated macrophages. If*ϕ*is the fraction of total macrophages that are in the activated state, (*ϕ*=*M*_{a}/*M*_{tot}, 1−*ϕ*=*M*/*M*_{tot}), then*J*+(*bϕ*−c(1−*ϕ*))*M*_{tot}≈0 at the threshold of chronic inflammation. We can explore various assumed values of 0≤*ϕ*≤1 and 1×10^{6}≤*M*_{tot}≤1×10^{8}to find corresponding values of*b*/*c*, and hence compute*b*using*c*≈0.1. For example, using realistic ranges such as*ϕ*=0.5,*M*_{tot}≈1×10^{7}corresponds to*b*≈0.09 d^{−1}, and this is the value used in our model. Since*M*_{tot}is physically limited to about 1×10^{8}cells ml^{−1}, reasonable anti-crowding terms are*e*_{1}=*e*_{2}≈1×10^{−8}.

With estimates of parameters as above, the required existence of the healthy rest state for both strains implies that(5.1)(using Balb/c parameters, or *k*/(*ℓ*−*k*)>10 using NOD parameters). However, a state of chronic inflammation for NOD mice only exists when(5.2)But these two conditions are contradictory, and cannot both be satisfied: either all mice invariably develop chronic inflammation (if for Balb/c or <10 for NOD), or an apoptotic wave never triggers chronic inflammation (in the case ). Furthermore, for NOD parameters, our estimates imply , violating condition (4.4*a*) by an order of magnitude. This points to our conclusion that NOD mice could *never* become inflamed in a biological range of parameters of the Copenhagen model, regardless of other strain-specific differences, e.g. in their values of *c*, *k*, or . Once quantified, it becomes evident that the disparities in macrophage function alone, in the context of the Copenhagen model, cannot explain the early chronic inflammation in NOD mice. We therefore ask what else could be missing, and consider the role of necrosis.

## 6. Necrosis as the inflammatory effect

In view of the above, we consider a revised model wherein damage results from interaction between activated macrophages and cells in secondary necrosis accumulating from insufficient clearance of apoptotic cells. Clearance by macrophages of cells undergoing apoptosis is generally associated with anti-inflammatory responses, while macrophage responses to necrotic cells, including secondarily necrotic cells derived from uncleared apoptotic cells, are perceived as pro-inflammatory (Gregory & Devitt 2004). We therefore assume that, upon encountering necrotic cells (*B*_{n}), activated macrophages (*M*_{a}) produce (at rate *α*) some harmful factor(s) (such as the cytokine IL1-β, see Stoffels *et al*. (2004), iNOS, or other chemicals). We represent such toxic factors collectively as a single variable (*C*), and assume a linear rate of removal (*δ*). Based on fig. 1 in Eizirik & Mandrup-Poulsen (2001), we further assume that the rate of cytokine-induced apoptosis of β-cells is a Michaelis–Menten saturated function of *C* with maximal rate *A*_{max} and half-max cytokine concentration *k*_{c}. (Saturation is essential to avoid unrealistically high damage during the neonatal apoptotic wave, when *C* could be elevated by 1–2 orders of magnitude.) To avoid introducing numerous arbitrary new parameters, we here assume that removal of necrotic cells by phagocytosis is identical to the removal of apoptotic cells (rates *f*_{1}, *f*_{2}). This extended model is shown in figure 2 and given in equations (6.1*a*)–(6.1*e*) (see table 3 for parameter values and definitions):(6.1a)(6.1b)(6.1c)(6.1d)(6.1e)We note that making a QSS assumption on the toxic factor (*C*′(*t*)≈0) leads to *C*_{QSS}≈(*α*/*δ*)*B*_{n}*M*_{a}, and reduces the equation for to(6.2)where . Effectively, this means that the previous damage rate per macrophage, , has been replaced by the term (*A*_{max}*B*_{n})/(*k*_{b}+*B*_{n}*M*_{a}), with two new parameters, A_{max} and *k*_{b}. Further assuming a QSS on *B*_{n} reduces the model to a system consisting of the original equations for *M*, *M*_{a}, together with the following equation for *B*_{a}:(6.3)The significant difference here is that the second term, corresponding to damage-induced apoptosis of β-cells, is very small when both *B*_{a} and *M*_{a} are small. Hence, for a low level of inflammation in the pancreas, the positive feedback leading to self-destruction of β-cells is negligible. (This addresses an unrealistic aspect of the basic Copenhagen model.) Consequently, in the model with necrosis effects, the healthy rest state is always stable. Further, the ‘effective damage function’ shown in equation (6.3) saturates with respect to both *B*_{a} and *M*_{a}, with maximal rate of apoptosis induction *A*_{max}. This means that the highest damage rate *per* *M*_{a} lies at intermediate inflammation levels, so β-cell destruction is bounded. (Otherwise, if this were not the case, a very rapid and complete destruction of the β-cell mass could occur even before the adaptive immune system gets involved. But this is never observed in NOD mice, and therefore cannot be realistic.)

Ideally, if we could identify the mechanism of induced damage, and dose–response behaviour of apoptosis caused by cytokines or other cytotoxins, we could base the new parameters *α*, *δ*, *A*_{max}, and *k*_{c} on such details, but this is premature at this stage of our work. Hence we use a ‘generalized’ chemical *C*, whose level at inflamed conditions far exceeds its half-max apoptosis-inducing level, *k*_{c}. To parametrize this extended model, we make the following estimates.

In the steady state corresponding to chronic inflammation in NOD mice, we take as reasonable cell densities, and (*B*_{a}+*B*_{n})≈2×10^{5}cells ml^{−1}. This means that cells ml^{−1}, or about 2.5% of the β-cell mass enter apoptosis each day. At this steady state, production and removal of apoptotic β-cells must balance, so that *A*_{max}*C*/(*k*_{c}+*C*)≈1×10^{7}. Assuming that the maximum β-cell-mass destruction due to the native immune system is around 5% per day, means *A*_{max}≈2×10^{7} cells ml^{−1} d^{−1}. We take the cytokine concentration relative to an arbitrarily chosen level for half-max apoptosis, *k*_{c}=1 nM. Only the relative levels of maximal *C* and *k*_{c} are important to the model, with the proviso that the Michaelian apoptosis term should be saturated during the peak of the apoptotic wave, and half-saturated during the chronic phase. The ratio *α*/*δ* determines cytokine accumulation. Assuming a short cytokine decay time of 1 h, i.e. *δ*≈25 d^{−1}, the cytokine secretion rate should be roughly *α*≈5×10^{−9} nM cell^{−2} d^{−1}. Then, in NOD mice during the chronic stage, , i.e. at the required level for inducing a half-maximal apoptosis rate.

We simulated equations (6.1*a*)–(6.1*e*) with parameter values indicated in tables 1 and 3, for both NOD and Balb/c parameters. (Using (6.3) instead of (6.1*c*)–(6.1*e*) leads to quantitatively indistinguishable dynamics, i.e. the QSS approximation is generally valid.) The two strains were assumed to differ *only* in the values of macrophage activation and phagocytosis rates, *g*, *f*_{1}, *f*_{2}, as indicated in table 1. Initial conditions were chosen to simulate a healthy rest state, i.e. *M*=4.77×10^{5}, *M*_{a}=0, *C*=0 and *B*_{n}=0. We incorporate the neonatal apoptotic wave *W*(*t*)=4×10^{7} exp(−((*t*−9)/3)^{2}) cells ml^{−1} d^{−1}, which means that the wave peaks after 9 days, when 4×10^{7} cells ml^{−1} d^{−1} enter apoptosis. Log-density time-plots are shown in figure 3. In NOD mice, the wave of apoptotic β-cells leads to a huge elevation of activated macrophages, up to 5×10^{6} cells ml^{−1}, and significant increase in resting macrophages, to roughly 1.1×10^{6}, after a phase of rapid activation of most resting macrophages. The level of apoptotic β-cells declines right after the normal remodelling stage, but stays elevated, at roughly 2×10^{5} cells ml^{−1}, due to inflammation-induced damage. In contrast, Balb/c mouse macrophages are able to clear the apoptotic wave rapidly enough that chronic inflammation is prevented. In Balb/c mice, activated macrophages, apoptotic β-cells and other byproducts (*C*, *B*_{n}) decrease to zero and resting macrophages are restored to their normal level in the tissue within about 10 days after the apoptotic wave ends.

## 7. Discussion

In this paper, we explored the hypothesis that T1D is initiated from what should be normal developmental β-cell turnover (Trudeau *et al*. 2000; Mathis *et al*. 2001). For a parsimonious explanation (T1D often found in NOD mice, never in Balb/c mice), we asked whether macrophage defects *alone*, could suffice, and whether the wave of apoptosis can trigger chronic inflammation. This phase must correspond to limited death of β-cells (a requirement of our model), which sets the stage for eventual T-cell priming associated with full-blown disease (not described here). We started with the qualitative ‘Copenhagen’ model by Freiesleben De Blasio *et al*. (1999), but found that it cannot be quantitatively correct within biologically reasonable parameter regimes even if other strain-specific variations are included (e.g. differences in macrophage-induced rate of damage to β-cells).

We extended the model in a realistic way, by assuming that apoptotic β-cells are cleared ‘silently’, whereas necrotic β-cells trigger secretion of damaging cytokines by activated macrophages (Stoffels *et al*. 2004). Since secondary necrosis is a post-apoptotic stage, inefficient clearance of apoptotic β-cells by macrophages gives rise to prolonged chronic inflammation for NOD, but not Balb/c mouse parameters, as shown in figure 3. The extended model has a ‘generic cytokine’ (e.g. IL-1*β*) that induces apoptosis with saturating kinetics: a saturation assumption is essential to avoid unrealistic rapid destruction of the β-cell mass. This is an important feature of any model for *early* events, as T-cells are known to be required in T1D pathogenesis.

Our parameter estimates are taken from previous experiments reported in Marée *et al*. (2005), and reasonable ranges for cell populations and turnover rates from other papers. (Our value for macrophage turnover rate *c*≈0.1 d^{−1} is based on work by Van Furth & Diesselhoff-den Dulk (1984), and differs from the value cited by Wigginton & Kirschner (2001) by an order of magnitude.) Our extended model generates appropriate dynamics in a wide range of parameter values. For parameter values such as *c* and *k*, whose estimated values vary widely, appropriate dynamics are found over the whole range of possible values. Extended bifurcation analysis (not shown here) points to other interesting dynamics, including limit cycle dynamics within certain parameter ranges. This will be the subject of a mathematical treatment elsewhere.

Clearly, many other factors could be at play in the differences between NOD and non-diabetic mice. Stoffels *et al*. (2004) suggest that NOD mice fail to downregulate apoptosis-inducing cytokines such as TNFα and IL-1β in response to apoptotic cells. Thus, rates of damage or secretion (, or *α*) could differ in the two strains as well. We have already argued that this would not rescue the original Copenhagen model; however, such effects would accentuate the ability of the revised model to account for the tendency of NOD, but not Balb/c to get the disease.

Future efforts should include experiments to check and/or refine parameter estimates, identify specific cytokine levels *in vivo* and their effects on β-cell apoptosis, and, eventually, include participation of T-cells in a comprehensive model for T1D. Moreover, the anti-inflammatory response to clearance of apoptotic cells versus pro-inflammatory response to clearance of necrotic cells should be better quantified, since this difference forms the basis of our explanation for why NOD, but not Balb/c, mice potentially develop diabetes. It can be studied either by *in vitro* engulfment experiments, or by inducing different levels of necrosis *in vivo* and monitoring cytokine levels. Such quantitative studies can also help to determine why inflammation is generally restricted to the pancreas. For example, in the developing rat kidney, an apoptotic wave has also been observed (Coles *et al*. 1993), but that wave does not trigger chronic inflammation. The neonatal apoptotic wave of β-cells has now been observed in several species, including mice and rats (O'Brien *et al*. 2002*a*), and there are strong indications for its presence in humans as well (Kassem *et al*. 2000). Moreover, impaired phagocytosis by macrophages of diabetic patients has been established (Katz *et al*. 1983; Geerlings & Hoepelman 1999). Here we have shown that incomplete clearance of an apoptotic wave can lead to chronic inflammation that sets the stage for the development of T1D.

## Editors' note

Please see also related communications in this focussed issue by Bernard *et al*. (2006) and Leem *et al*. (2006).

## Acknowledgments

We gratefully acknowledge support by the Mathematics of Information Technology and Complex Systems (MITACS). A.M. is supported by the Research Council for Earth and Life Sciences (ALW) with financial aid from the Netherlands Organization for Scientific Research (NWO). R.K.'s graduate fellowship is sponsored by the Natural Sciences and Engineering Research Council, NSERC (Canada). D.F. is also supported by the Canadian Institutes for Health Research. We thank Marek Łabecki for initial steps towards this modelling effort, and Cheryl Dyck and Mitsuhiro Komba for helpful discussions.

## Footnotes

One contribution of 13 to a Theme Issue ‘Biomathematical modelling I’.

- © 2006 The Royal Society