## Abstract

Macrophages play an important role in clearing apoptotic debris from tissue. Defective or reduced clearance, seen, for instance, in non-obese diabetic (NOD) mice, has been correlated with initiation of autoimmune (Type 1) diabetes (T1D) (O'Brien BA, Huang Y, Geng X, Dutz JP, Finegood DT. *Diabetes* 51: 2481–2488, 2002). To validate such a link, it is essential to quantify the reduced clearance (for example, by comparison to BALB/c control mice) and to determine which elements of that clearance are impaired. Recently, we fit data for the time course of in vitro macrophage feeding experiments to basic models of macrophage clearance dynamics, thus quantifying kinetics of uptake and digestion of apoptotic cells in both mouse strains (Marée AFM, Komba M, Dyck C, Łabeçki M, Finegood DT, Edelstein-Keshet L. *J Theor Biol* 233: 533–551, 2005). In the cycle of modeling and experimental investigation, we identified the importance of *1*) measuring short-, intermediate-, and long-time data (to increase the accuracy of parameter fits), and *2*) designing experiments with distinct observable regimes, including engulfment-only and digestion-only phases. Here, we report on new results from experiments so designed. In comparing macrophages from the two strains, we find that NOD macrophage engulfment of apoptotic cells is 5.5 times slower than BALB/c controls. Significantly, our new data demonstrate that digestion is at least two times slower in NOD, in contrast with previous conclusions. Moreover, new data enable us to detect an acceleration in engulfment (after the first engulfment) in both strains, but much smaller in NOD macrophages.

- NOD mice
- autoimmune diabetes
- macrophage phagocytosis
- engulfment rate
- digestion rate

a variety of innate and adaptive immune system components are implicated in the destruction of β-cells in the pancreatic islets of Langerhans, causing autoimmune diabetes [Type 1 diabetes (T1D)]. Self-reactive T cells are known to be the main effectors of β-cell killing (10, 19). However, controversy exists as to the events that initiate the activation of β-cell-specific T cells (9, 35). Macrophages are among the earliest cells to infiltrate the islets. In the nonobese diabetic (NOD) mouse, a common animal model for T1D, this infiltration occurs as early as 2 wk of age and coincides with a normal wave of apoptosis of β-cells (40). This is followed by infiltration of dendritic cells and CD4^{+} and CD8^{+} T cells at 6–10 wk (8). Female NOD mice become diabetic by 12–15 wk (19) when 95% of their β-cells have been destroyed (23).

Relative to BALB/c (controls), macrophages of NOD mice are identifiably different (16, 39), exhibiting abnormal cytokine secretion (2), reduced activation of other immune cells (17, 20), and reduced phagocytosis (35). Based on these observations, it has been suggested that defective macrophage clearance of the neonatal wave of β-cell apoptosis leads to self-presentation and, subsequently, to T-cell activation (27, 40). In earlier mathematical models, it has been shown that reduced macrophage clearance could be a major factor in T1D initiation (12, 25). These observations have motivated an ongoing comparative study of macrophage behavior in NOD vs. BALB/c mice.

We recently showed that interstrain differences can be inferred by fitting time-course data of in vitro macrophage feeding experiments to mathematical models of macrophage dynamics (25). Our previous study demonstrated the feasibility of determining multiple kinetic parameters from such data. Moreover, it identified the importance of measuring data on multiple time scales and at different regimes: one of the improvements described here is inclusion of a “digestion only” regime to directly measure the rate that macrophages process apoptotic bodies after engulfment. Our main goal in this paper is to report on results of these redesigned experiments and to show how these facilitate more accurate determination of clearance kinetics and, hence, better inter-strain comparison.

Before comparing between strains, we ask basic general questions about macrophage clearance. *Question 1*: After the first engulfment, do subsequent engulfments take place at higher rate? If so, we refer to the acceleration in engulfment as an “activation” step.^{1} *Question 1a*: If activation is observed, is that activation reversible or irreversible on the time scale of several hours? *Question 2*: Is there evidence for a maximal capacity of engulfed apoptotic bodies in a macrophage? *Question 2a*: If so, what is that maximal number? *Question 3*: Are the engulfed apoptotic bodies digested simultaneously, or does digestion saturate?

In *question 1a*, we distinguish between the case that a macrophage reverts to a quiescent state (reversible activation) after fully digesting a meal or remains in an elevated engulfment mode (irreversible activation) over the time scale of a few hours. In *question 3*, we distinguish between *1*) digestion that can handle many apoptotic bodies in parallel, *2*) a saturated digestion machinery that can only break down one body at a time, and *3*) a machinery that saturates only at higher numbers of apoptotic bodies. (Note that we refer to apoptotic cells preengulfment and to apoptotic bodies while still visible inside a macrophage.) By comparing the quality of fit of data to several model variants, we distinguish between competing models and answer (*questions 1*–*3*). Having identified the optimal model (and hence the most likely macrophage dynamics), we assess differences between the control (BALB/c) and diabetes-prone (NOD) mice.

## MATERIALS AND METHODS

### Animals

Five-week-old female BALB/c and NOD/Lt mice (purchased from Taconic, Germantown, NY) had free access to standard laboratory chow and water. Before collection of peritoneal cells, mice were killed by CO_{2} asphyxiation. Animals were handled in accordance with the regulations of the Canadian Council on Animal Care, and all experimental manipulations of the animals were approved by the Animal Care Committee at Simon Fraser University.

### Isolation and Culture of Macrophages

Nonelicited macrophages were collected from mice by peritoneal lavage with cold complete medium (RPMI 1640; 100 U/ml penicillin, 100 μg/ml streptomycin, and 10% vol/vol heat-inactivated FCS, all from Life Technologies, Burlington, Ontario, Canada). Macrophages from multiple mice of a given strain were pooled together, and aliquots were seeded into eight-well chamber slides (Nalge Nunc). Cells were incubated at 37°C/5% CO_{2} for 2 h, allowing macrophages to adhere and spread on the slides. Unbound cells were washed three times with PBS, then cultured with apoptotic cells. The resulting number of adherent macrophages was ∼1 × 10^{5} cells per well.

### Induction of Apoptosis

A suspension of thymocytes harvested from BALB/c mice was irradiated by exposure to 254-nm ultraviolet C light for 20 min, and then it was cultured in serum-free complete medium for 2 h at 37°C/5% CO_{2} to let apoptosis progress. FCS (10% vol/vol) was then added to the medium.

### Phagocytosis and Digestion Assays

Apoptotic cells were cultured with macrophages at a ratio of 5:1 apoptotic cells-to-macrophages at 37°C/5% CO_{2} for variable periods of time. A total of six wells were used for each strain of mice for each experiment. Each well was assigned a specific time point: 20 min, 1 h, 2 h, 3 h, 4 h, and 5 h. After 20 min, 1 h, and 2 h, cells from one well per strain of mouse were washed three times with PBS to remove apoptotic thymocytes that had not been phagocytosed. For the “engulfment” experiment, cells were then fixed in 10% neutral buffered formalin (Sigma, Oakville, Ontario, Canada) and stained with Mayer's haematoxylin and eosin (Sigma) for identification of apoptotic thymocytes and macrophages.

For the “digestion-only” phase, at the end of 2 h of co-culture of macrophages and apoptotic cells, the remaining three wells per strain were washed with PBS to remove residual (unengulfed) apoptotic cells. Complete medium was then added to the wells, and adherent cells were further cultured for 1, 2, or 3 h (for the time endpoints at 3, 4, and 5 h, respectively, on ⇓Fig. 2). At the end of a culture period, cells were washed, fixed, and stained as above. Cells were observed with light microscopy at 1,000-fold magnification under oil immersion.

Phagocytosis was evaluated by counting and classifying ∼1,000 macrophages per well. Altogether, 60 wells (30 per strain) were analyzed, with a total of 58,646 macrophages. Each macrophage was assigned to a class (*i* = 0, 1, 2,… *N*) according to the number, *i*, of internalized apoptotic bodies it contained. Generally, up to *N* = 7 internalized apoptotic bodies could be observed, but sporadically even larger numbers, up to 11, were seen. Only apoptotic bodies clearly visible within the perimeter of the macrophage were counted. Results were initially expressed as the fraction of macrophages within a certain class in each well. These numbers were subsequently converted into macrophage densities (in units of cells/ml) by multiplying the fraction by the number of adherent macrophages per well (1 × 10^{5}) and dividing it by the well volume (0.5 ml).

### Modeling Techniques

Our modeling techniques are based on Marée et al. (25). Definitions of all variables and parameters are given in Tables 1 and 2, and a full composite of the models is shown in Fig. 1. These include phagocytosis with and without an activation step (which is either reversible or irreversible), three variants of digestion (parallel, serial, and saturating), and two distinct maximal number of engulfed bodies per macrophage (*N* = 7 or 12). In each case, *M*_{i}(*t*) represents the number of macrophages that contain *i* visible apoptotic bodies at time *t* (hours), with *i* = 0, 1, 2, etc. up to *N*, and *A* is the remaining number of unengulfed apoptotic cells. Macrophages “transit” between classes by engulfing apoptotic cells (arrows pointing *right* in Fig. 1), using up *A*, or by digesting an internalized cell (arrows to *left*). On the time scale of the experiment (5 h), the macrophage population is approximately constant (2 × 10^{5} cells/ml).

We simplify the macrophage clearance (targeting, binding, engulfing, and digesting apoptotic or necrotic cells) as essentially two processes, each with its own inherent time scale, engulfment [represented by parameter *k*_{e} (ml·cell^{−1}·h^{−1})], and digestion [represented by parameter *k*_{d} (h^{−1})]. Experimental evidence (11, 22, 35) suggests that, after the first engulfment, subsequent engulfment are either faster (22, 35) or slower (11) than the initial rate. To allow for either case, we use *k*_{a} (ml·cell^{−1}·h^{−1}) to represent the initial engulfment rate and *k*_{e} for the later rates. (Because we find that *k*_{a} < *k*_{e}, we later refer to *k*_{a} as the “activation rate” and the basal engulfment state before activation as the “rest state.”) We considered two cases: *1*) macrophages stay activated on the experimental time scale vs. *2*) return to rest state after fully digesting a meal [see, e.g., Licht et al. (22) who observe higher engulfment rates even after 20 h].

We considered several possible maximal capacities of apoptotic bodies per macrophage, *N*. In one case, *N* was set to the maximum observed experimentally; we also tested a much larger value, and a third case where such a limit does not exist. Since details of the digestion process are not known, we considered the following possibilities: *1*) a saturating Michaelian rate of the form *k*_{d}*i*/(c + *i*), where *i* is the load (i.e., number of apoptotic bodies to be digested), and c is the load at half-saturation; *2*) a linear regime, with rate (*k*_{d}/c)*i* = *k*_{d}*i* (which approximates *regime 1*) if saturation occurs only at a very high load (this is termed “parallel digestion” since the time to digest any load *i* is the same); *3*) serial digestion, i.e., the opposite regime where the process is fully saturated, and the digestion rate is *k*_{d}.

To avoid combinatorial proliferation of possible models, we here discuss a subset of models chosen to explore the main features of interest (Fig. 1).

#### Basic model.

The engulfment rates are all identical (*k*_{e0} = *k*_{e}), and digestion is serial with identical constant rates (*k*′_{d} = *k*_{d2} = *k*_{d3} = *k*_{d}) (see Fig. 1*A*).

#### Variants I and II: reversible/irreversible activation.

The first and subsequent engulfments occur at distinct rates (*k*_{e0} = *k*_{a} ≠ *k*_{e}). In variant I, we identify the first step as a reversible activation step, which means that digesting the last apoptotic body brings the macrophage back to the rest state (Fig. 1*A*). In variant II, we included a separate class, *M*_{0a}, of activated macrophages without apoptotic bodies. Only the resting macrophages, denoted *M*_{0r} in this variant, have a distinct engulfment rate *k*_{e0} = *k*_{a} ≠ *k*_{e} (Fig. 1*B*). The classes *M*_{0a},*M*_{0r} cannot be distinguished experimentally but can be determined by fitting the data to the model.

#### Variants III and IV: parallel/serial digestion.

In variant III, we modified variant II by assuming parallel digestion, i.e., *k*_{d1} = *k*′_{d}, *k*_{d2} = 2*k*′_{d}, *k*_{d3} = 3*k*′_{d}, *k*_{di} = *ik*′_{d}. In variant IV, we assumed Michaelian digestion rates, *k*_{d1} = *k*_{d}/(c + 1), *k*_{d2} = 2*k*_{d}/(c + 2), *k*_{d3} = 3*k*_{d}/(c + 3), *k*_{di} = *ik*_{d}(c + *i*), where *i* is the number of apoptotic bodies to digest and c is the number at which digestion is at half of its maximal rate.

### Parameter Fitting and Model Comparison

As in Marée et al. (25), we fit parameters in two ways: first, by using simple features of the data, and second, by full fitting using the freeware biochemical simulator program Gepasi and its Levenberg-Marquardt optimization routine (29–31).

Because cells have to be fixed before counting, the closest we can get to a true time series is to run experiments concurrently, using a single aliquot of macrophages for all wells. We call such a time series “a single replicate.” In a concurrent fit, we fit to all data points of all replicates simultaneously, although they are obtained from different aliquots of macrophages. Here, a confidence interval for parameter values can be calculated based on the drop in the goodness-of-fit when each individual parameter is changed. (Alternately, fitting each replicate separately preserves the distinction between aliquots.) With replicate fitting, one can calculate the standard deviation in the default way, because the outcome of the fittings can be considered to be five independent measurements of the same kinetic constant. By comparing the quality-of-fit between concurrent and replicate fitting, we can also determine whether there are significant differences between replicates (which would warrant critical evaluation of the experimental setup). A full discussion of the advantages and disadvantages of each method is provided by Marée et al. (25).

To select the model that best fits the observed data, we used the Akaike's Information Criterion (AIC) (1, 32). The criterion selects the model that has greatest accuracy with fewest free parameters^{2} [see appendix and fuller discussion in Marée et al. (25)]. The Akaike weights, *w*_{i} (Table 5), represent the relative probability of each candidate out of a set of models being the best representation of the dynamics underlying the data. We also used this criterion to compare concurrent vs. replicate data fits.

## RESULTS

We describe parameter values from ballpark estimates and from full data fits. Our procedure is explained below, and the results are summarized in Tables 3 and 4. The last column in these tables is the predicted fold difference in the parameter values between NOD and BALB/c mice and the 95% likelihood interval calculated using the default Student's *t*-test to compare two samples with different standard deviations. An asterisk indicates *P* < 0.05.

### Ballpark Parameter Estimates

Our models use features of experimentally observed dynamics of macrophage classes (Table 1) to compute rough estimates of parameters for individual macrophage function (Table 3). The relevant formulae are based on simple features of the data such as slopes ([s0, s1, s2]), ratios ([r1]), and other descriptors ([p1], [i1]), as summarized in Table 2 and cited by tags in square brackets here. [For full derivations, see Marée et al. (25).] For accurate parameter fits, observations are needed at early (*t*_{1} = 20 min), and intermediate (*t*_{Q} = 1 and 2 h) time points, as well as in the digestion phase (*t*_{d} = 2, 3, 4, and 5 h after washing out remaining apoptotic cells). The quality of the estimates improves if *A* is also measured at each time point (not possible in our experiments).

#### Early time points.

In the beginning of the experiment, all macrophages are in class 0 (*M*_{0} = *M*). The dominant process, engulfment, is then related to the rate of decrease of *M*_{0} and/or rate of increase of *M*_{1}, i.e., to the initial slope of the data curve. An approximation of this derivative using a finite difference of data values (at discrete times, e.g., *Eq. 3*) is reasonable, provided that time points are close together. This early time data is thus used to estimate the engulfment rates [*k*_{e} in the basic model ([e1,e2]), and *k*_{a} in variants I and II ([a1])]. In the basic model, all subsequent engulfments are identical. In activation-variant models, the rates of (subsequent) engulfments can be obtained by rates of change of higher classes at early time points ([e3]). Time points too far apart underestimate these parameters.

#### Distribution of macrophage classes and parameter ratios.

After about 1 h of feeding, the number of apoptotic cells and ability of macrophages to engulf and digest equilibrate. This approximate equilibrium is termed a quasi-steady state (QSS). At this time (*t* = *t*_{Q}), ratios of the sizes of successive classes are roughly constant ([r1]), that is *M*_{1}/*M*_{0} ≈ *M*_{2}/*M*_{1} ≡ *R* = constant. Furthermore, in the basic model, it can be shown that this ratio is related to a ratio of model parameters: (1) If *A* is observed experimentally, this ratio can be used to calculate *k*_{d} based on the previously obtained *k*_{e} value ([d1]).

#### Estimates for the digestion rate.

Studies of phagocytosis often cite the percent phagocytosis Φ (the percent of macrophages that have engulfed one or more apoptotic cells) and the phagocytic index *I*_{Φ} (the number of visible apoptotic bodies per 100 macrophages) (33, 35). These indexes ([p1, i1]) are helpful in obtaining an estimate of the digestion rate, and they are also related to the ratio *R* ([r3, r4]).

In the digestion phase, apoptotic cells have been washed out (*A* = 0), and no further engulfments occur, so each macrophage gradually digests its meal and “transits” to a lower class. The number of macrophages in a certain class could still increase due to a large influx of macrophages from a higher class, but the total number of apoptotic bodies within all macrophages should continuously decrease. In the case of parallel digestion, the decrease in total number of apoptotic bodies should be exponential. Thus the phagocytic index *I*_{Φ} (which is the relative total number of apoptotic bodies) should also decrease exponentially (*Eq. 4* in appendix). Hence, it follows that (2) Here ln( ) is the natural logarithm, *Ĩ*_{Φ} is the phagocytic index at the start of the digestion phase (*t* = *t*_{d0} = 2 h), and *T* is elapsed time (*T* = *t* − *t*_{d0}) (see formula [d4]).

A number of other formulae ([d2, d3, e4, e5]) link various ratios and parameters. All have been derived for our specific set of models, but these can be generalized to other situations.

### Model Comparison by Full Data Fits

Because the above estimates are fairly rough, for model comparison we use full data fits. We compare all model variants (and concurrent vs. replicate data fits) in Table 5 to select the best model. *Rows 1*–*12* compare models for BALB/c data, and *rows 13*–*24* compare models for NOD data. The last 12 rows compare digestion variant models using only the digestion data.

#### Macrophage activation.

We find that models with an activation step are more likely than the basic model and that consecutive engulfments are indeed faster than the initial engulfment. [*k*_{e} is significantly different and always higher than *k*_{a}, contrary to reports in Erwig et al. (11)]. Reversible and irreversible activation are found to be equally likely (variants I and II; Table 5). The fact that activation is an important feature of the behavior in both strains is a new result that differs from our previous conclusions.

#### Macrophage capacity.

In variant II, comparing *N* = 7 to *N* = 12 results in no difference in fit. If there were an underlying maximum capacity, the model with *N* = 12 that predicts higher levels of apoptotic bodies would have a lower goodness-of-fit and be less likely. Since the same dynamics (and therefore the same likelihood) are predicted for *N* = 12 and higher *N*, a phagocytosis ceiling did not play a role in our experiments. (If such a ceiling exists, a much larger meal of apoptotic cells would be needed to detect it; this would fill higher *M*_{i} classes.) Our data indicate that the observed maximum number of apoptotic bodies per macrophage (i.e., 7) is not due to limitations on further engulfment.

#### Serial, parallel, or saturating digestion.

Serial digestion is most probable for both strains in each case. Although the SSE is lower for the saturating model (variant IV), the improvement of fit is too small to justify an extra parameter, so serial digestion most parsimoniously explains the data. However, for NOD macrophages, parallel and saturating digestion (variants III and IV) still have significant probability of best describing the data. Because NOD mice have strongly reduced phagocytic capabilities, only a small fraction of their macrophages are in higher classes. Thus it becomes harder to compare digestion rates for different classes, as is needed to discriminate between the digestion model variants. Altogether, our model comparisons indicate that digestion is best described as a serial process in both mouse strains.

#### Replicate or concurrent fitting.

The likelihood of the replicate models is essentially zero. That is, we find no indication for variation in experimental setup or performance between replicates. There is effectively no replicate-to-replicate variation in the kinetic constants. This observation allows us to extrapolate from the outcome of this study (i.e., the estimated parameter values) to more general quantitative questions about the dynamics underlying the development of diabetes in NOD mice [see Marée et al. (25)].

### Parameter Values Obtained

The results of full fitting of models to concurrent data using the optimization feature of the software Gepasi are shown in Table 4. Clear differences emerge between the parameters for control (BALB/c) vs. diabetes-prone (NOD) mice, revealing details not apparent from the ballpark estimates of Table 3. As discussed in model comparisons, the parameter values given by model variants I and II should best estimate macrophage clearance dynamics, with variant II (*N* = 7) having the highest likelihood for both strains. Below, we compare across models and strains. Values are given for variant II (*N* = 7) unless otherwise indicated.

#### Engulfment rates.

The initial engulfment rates *k*_{a} are 2.59 ± 0.03 times greater for BALB/c than for NOD macrophages, and subsequent engulfment rates (*k*_{e}) differ more dramatically (by 5.54 ± 0.12; see Table 4). Model variants I and II (*N* = 12) give highly comparable parameter estimates. This confirms that NOD macrophages are strongly impaired in engulfment. The rough estimates for *k*_{a} (variant I; Table 3) are both qualitatively and quantitatively close to the full data fit values. Qualitatively, the ratios of BALB/c to NOD *k*_{e} values are similar for both fitting methods, but absolute levels are lower in ballpark estimates than in full data fits. The engulfment rate estimated by the full fit of the basic model is between the initial and subsequent rate; the ballpark estimate simply gives the initial rate. In summary, the ballpark estimates give a reasonable representation of the engulfment dynamics.

#### Activation.

For BALB/c macrophages, Table 4 reveals a clear difference between the rate of engulfment of the first apoptotic cell *k*_{e0} = *k*_{a} = 5.4 × 10^{−7} ml·cell^{−1}·h^{−1}, and the subsequent, much faster engulfment rate ke = 23.6 × 10^{−7} ml·cell^{−1}·h^{−1} (more than five times higher). There is also a small but significant activation step in NOD, from *k*_{a} = 2.1 × 10^{−7} ml·cell^{−1}·h^{−1} to *k*_{e0} = *k*_{a} = 4.3 × 10^{−7} ml·cell^{−1}·h^{−1}. In the ballpark estimates, increases in the engulfment rates are strongly underestimated and hardly reveal any activation in the NOD case. (Full data fitting becomes essential for such fine detail.) Importantly, in contrast with Marée et al. (25), the fact that NOD macrophages also undergo activation could only be observed using both the improved dataset and full parameter fitting.

#### Digestion rates.

The rate of digestion of the internalized apoptotic cells are 2.41 ± 0.05 as large in BALB/c as in NOD. [For model variant I, a smaller fold-difference (1.97 ± 0.05) is observed.] This means that digestion is impaired in NOD macrophages, but not as dramatically as engulfment. The number of apoptotic bodies at half-saturation (variant IV) is c ≈ 1 for NOD and lower for BALB/c: in both strains, a single apoptotic body nearly saturates digestion. This explains why the experimental data are much better described by serial digestion (the immediate saturation extreme) than by parallel digestion (the no saturation extreme). The ballpark digestion rates are quite similar but clearly overestimated. Only the parallel digestion estimate (variant III), based on the digestion-only phase, is similar to the full data fit. Together, this shows that for good ballpark estimates of *k*_{d}, it is essential to include a digestion phase in the experiments (unless the number of remaining apoptotic cells, *A*, is measured during the experiment). Otherwise, the accuracy of *k*_{d} is quite limited as a consequence of being estimated indirectly. The more accurate data for the digestion phase and full data fits revealed a previously undetected interstrain difference.

### Displayed Dynamics

Figures 2 and 3 show the dynamics of the macrophage engulfment models fit to data. We show the most likely model (variant II, *N* = 7) and compare it with the basic model (without activation) and variant IV (which gave the best but not the most parsimonious fit). The performance of BALB/c macrophages (*left*) is visibly more vigorous than NOD (*right*). A comparison of left panels shows that a model with an activation step fits BALB/c data significantly better than the basic model (compare black curves in each case). In the NOD fits, the disparity between models with and without activation is not as strong, but the basic model is also not as accurate as the models with an activation step (compare green curves). These results contrast with Marée et al. (25), where we found no significant improvement in the likelihood of a model with activation for NOD macrophages. Visually, no improvement of the fit can be observed for variant IV, so the usage of an extra parameter is not justified.

Figure 3 shows the phagocytosis index (*I*_{Φ}) and percent phagocytosis (Φ) time courses. These curves were not fit separately and give an independent indication of relative performance of models against data. Squares and triangles are Φ/*I*_{Φ} values computed from each data point, and curves are computed from model predictions using the best-fit parameters in Table 4. The top left panel confirms that the basic model does not account for the behavior of BALB/c macrophages, since the black curve representing *I*_{Φ} is significantly lower than the data points. In NOD mice (right panels), the distinction is not as great.

We computed the remaining unengulfed apoptotic cells, *A*, that could not be directly observed experimentally (simulations not shown). We found a huge difference systematically predicted between BALB/c and NOD mice: in 2 h, BALB/c macrophages clear 30% of the available apoptotic cells compared with <10% for NOD mice.

## DISCUSSION

### Major Conclusions and Insights

The results in this paper help to more accurately elucidate quantitative differences between macrophages in NOD mice prone to a form of autoimmune (Type 1) diabetes compared with control (BALB/c) mice. In agreement with previous work, we found that NOD mouse macrophages have a slow rate of engulfment of apoptotic cells compared with controls (25, 35). (Even after activation, engulfment rate in NODs is less than one-fifth of controls.) Like Marée et al. (25), but unlike older papers, we here quantified biological parameters that represent the actual rates of phagocytosis of apoptotic cells and the rate of digestion of those cells once they are engulfed. In contrast with Marée et al., a new major conclusion, based on the new refined experiments, is that activation is an important feature in both strains. An extension of Marée et al. and a second major conclusion is that digestion is serial (rather than saturating or parallel).

Our models informed experimental design twice (see Table 6 for an overview). First, models indicated that early time data is essential for accurate parameter estimates because intermediate time data only reveal quasi-steady-state ratios (such as *R*) rather than individual parameter values. Thus inaccurate engulfment rates engender poor estimates for digestion rates, since these parameters are linked. Furthermore, without early data, the distinction between models with and without an activation step is blurred. Second, even with early data (25) but without observing remaining apoptotic cells directly, a “ridge” in the parameter landscape hindered optimization of parameter fits. We then redesigned experiments to include early, middle, and late time behavior as well as a separate digestion phase; this led to proper estimation of the digestion rates and, hence, all linked parameters. A third major conclusion of this work is that, unless *A* can be observed over the course of the experiment, it is important to obtain data from a digestion-only phase of the experiments.

The various ballpark estimates link population-level phenomena (distribution of macrophages into classes) to underlying cellular behavior (individual macrophage rates of activation, engulfment, and digestion of apoptotic cells). These estimates are less accurate than full fits but reveal some interesting features. For example, the percent phagocytosis, Φ, is simply proportional to the quantity we called *R*, which contains a ratio of engulfment and digestion rates. (The phagocytic index, *I*_{Φ}, is also related to *R*, but not as simply.) This suggests that simple ratios of Φ could be more easily interpreted than ratios of *I*_{Φ} (see [r2] vs. [r3] in Table 2) to compare values of *k*_{e}/*k*_{d} across strains.

From Tables 3 and 4, highly significant fold-differences between strains can be observed. (An asterisk in the last column indicates *P* < 0.05.) It is possible that the fitting method underestimates the standard deviation in the parameter values, since it does not take into account the possible covariance between parameters. (This could lead to less dramatic drops in quality-of-fit when parameters are varied concurrently.) However, even the rough estimates in Table 3 obtained from formulae in Table 2 demonstrate significant differences between NOD and BALB/c strains. These, and also comparisons of parameter estimates from different cohorts (data not shown), lead to much larger confidence intervals, but nevertheless the fold-differences remain highly significant.

Assessing the relative merits of rough (Tables 3) and full (Table 4) parameter fits, we found that full fits using optimization software achieve more accurate results, as expected, but that the two sets of results are in the same order of magnitude. Furthermore, the ratios of rate constants for control and NOD mice are roughly the same. This implies that ballpark estimates would suffice to compare two conditions or strains, whereas full fits improve the absolute values of the parameters.

For model comparison, the Akaike criterion has wide applicability, wherever the details of underlying mechanisms are not known in advance. Using this technique, with the combined dataset from this experiment and from the first 2 h of the experiments reported in Marée et al. (25), we found that the variants I and II (reversible and irreversible activation) provide the most probable candidates for both BALB/c and NOD. To properly distinguish between these two variants, one should perform experiments in which the macrophages are presented a second apoptotic meal after the first digestion phase. Similarly, to investigate whether there is truly a maximum number of apoptotic bodies that can be found within a single macrophage (as well as to determine its value), one should feed macrophages with a higher density of apoptotic cells than was used in the current experiments.

### Implications and Connections With Autoimmunity

A major difference between mouse strains is that, in NOD macrophages, the engulfment rate of apoptotic cells is abnormally low by a factor of between 2.6 and 5.5. A minor difference is in reduced degradation of apoptotic cells once they are engulfed. Here, we mention some of the implications of our results in the context of autoimmunity and pathogenesis of T1D (27, 40).

Aside from in vitro experiments, we have recently shown a systemic defect in NOD macrophages in vivo [in the peritoneum, the skin, and the thymus (34)]. Challenging the NOD mouse with ultraviolet irradiation to the skin results in an increased apoptotic cell load (relative to control strains) and increased antibodies to DNA in circulation. Macrophage defects of a similar type have also been inferred in other autoimmune animal models, e.g., mice prone to lupus (7, 37). Mice that lack a membrane tyrosine-kinase (c-mer) have defective macrophage phagocytosis and are prone to a form of autoimmunity that resembles lupus (7); a mutant protein in mice that masks phosphatidylserine on apoptotic cells inhibits the phagocytosis of those cells, which then leads to autoantibody production in the animals (3). All these experimental studies link autoimmunity and defective clearance of apoptotic cells.

A connection to human disease has also been made. Reduced clearance of apoptotic cells has been observed in human patients with systemic lupus erythematosus (18). Gaipl et al. (14) showed that apoptotic cells accumulate in germinal centers of the lymph nodes and in the skin of systemic lupus erythematosus patients after exposure to ultraviolet light. They showed that attraction signals for macrophages were lower in the sera of about one quarter of the patients and suggested that defective clearance of dying cells could explain the accumulation of nuclear autoantigens in those patients. [See also Gaipl et al. (13) for a recent review, but note the contrasting study by Reefman et al. (38) where clearance rates of apoptotic cells were not significantly different.]

There are a number of hypotheses about the causes underlying the link between phagocytic defects and autoimmunity. One suggestion is that residual apoptotic material undergoes secondary necrosis and results in inflammation due to release of endogenous chemical danger signals (15). Whether apoptotic cells also release such inflammatory signals is still controversial, but it has been shown recently that Jurkat cells release strong danger signals (high-mobility group box 1 protein) to the extracellular environment 30 h after induction of apoptosis with pharmacological agents such as staurosporine (4). Although not verified experimentally, if secondary necrotic β-cells release such danger signals, this could induce maturation of dendritic cells that then prime naive autoreactive T cells. Furthermore, in some circumstances, inflammatory stimuli appear to reduce macrophage phagocytic ability, for example, in rat alveolar macrophages exposed to prolonged low levels of interferon-γ (5); this could lead to a positive feedback loop, further decreasing clearance.

Macrophages are known to play a role in the digestion of DNA from apoptotic cells (28, 36). Our new finding in this paper, that there are differences in this process (although to a lesser degree than differences in engulfment) between diabetes-prone and diabetes-resistant mice, could also have implications in autoimmunity. It has been shown that mice whose macrophages lack DNase II in lysosomes exhibit an autoimmune phenotype, polyarthritis, similar to human rheumatoid arthritis. Kawane et al. (21) showed that, when macrophages cannot degrade DNA from apoptotic cells, they produce TNF-α, stimulating other cytokine production and culminating in chronic polyarthritis. Apoptotic cells (such as Jurkat cells) release DNA into their medium in vitro (6). DNA that escapes degradation or is released from macrophages could cause immune response to self-tissues and initiate autoimmunity. This could be yet another mechanism linking defects in macrophage digestion to the pathogenesis of autoimmune diabetes.

Using an extension and revision of the Copenhagen model (12) with parameter values determined in Marée et al. (25) and taking the inflammatory effect of necrotic β-cells into account, we have shown that inefficient clearance of apoptotic β-cells by macrophages can, by itself, give rise to prolonged chronic inflammation for NOD but not BALB/c mouse parameters (26). This chronic inflammation could then set the stage for the development of T1D. In a later modeling study (24), we showed that the difference in clearance rates of apoptotic debris (and thus of putative auto-antigen peptides) between mouse strains can account for the fluctuating dynamics of CD8^{+} T cells observed by Trudeau et al. (41) in the prediabetic phase in NOD mice but not in controls. When clearance rates are normal, an immune stimulus that primes T cells gets resolved to a baseline state after a short time, and no autoimmunity occurs. A reduction of the clearance rate by a factor of 2 or more, as found here for NOD mice, could lead to cycles in the prevalence of auto-reactive T cells in our model, with eventual elevated levels of effector cells, autoimmunity, and death of β-cells.

In view of the above, our results have implications that go beyond simple quantification of phagocytosis. Indeed, such links suggest that macrophage clearance plays an important role in preventing autoimmunity and that defects in that function could account for pathogenesis in T1D. Treating macrophage defects and enhancing the low engulfment rate of apoptotic material in diabetes-prone animals could be an additional strategy in retarding autoimmune diseases such as T1D.

## APPENDIX

### Details of Formulae

The initial slope of the data curve in the engulfment phase can be approximated by (3) A similar result holds for *S*_{0}, the magnitude of slope of the curve *M*_{0}.

In the digestion phase (assuming parallel digestion), the change in the phagocytic index, *I*_{Φ}, satisfies a simple linear first-order differential equation whose solution is (4) where *Ĩ*_{Φ} is the phagocytic index at the start of the digestion phase, and *T* = *t* − *t*_{d0}. Rearrangement of this formula leads to *Eq. 2*.

### Akaike Weights

The Akaike's information criterion (AIC) (1, 32) maximizes the probability that a candidate model generated the observed data by rewarding descriptive accuracy while penalizing an increase in the number of free parameters. We have used the corrected version, (5) where *n* is the number of observations, *K* is the number of fitted parameters, and SSE is the remaining sum of squared error after curve fitting. The model with the lowest AIC_{c} value is considered to be the best model with most parsimonious description of the data. The Akaike weight (*w*_{i}) represents the probability that model *i* is the best model, given the data and the set of candidate models: (6) where *M* is the total number of models, the minimum is taken over all candidate models, and where Δ_{i}(AIC_{c}) = (AIC_{c})_{i} − min(AIC_{c}). Relative probabilities of any two models are proportional to the weights computed for them, regardless of the total number of models being compared. The advantage of the AIC method is that it allows us to compare all of the models simultaneously as easily as a subset of these models.

## GRANTS

We gratefully acknowledge support by the Mathematics of Information Technology and Complex Systems and by the Juvenile Diabetes Research Foundation. A. F. M. Marée is supported by the Netherlands Organization for Scientific Research. L. Edelstein-Keshet is also supported by the Natural Sciences and Engineering Research Council (Canada). D. T. Finegood is also supported by the Canadian Institutes of Health Research.

## Footnotes

↵* A. F. M. Marée and M. Komba contributed equally to this work.

↵

^{1}Our term “activation” should not be confused with similar term applied to transition of macrophages to antigen-presenting cells or other changes in response to inflammatory stimuli.↵

^{2}By introducing more free parameters, it becomes easier to fit the data more closely, and therefore, AIC penalizes each additional parameter.The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2008 the American Physiological Society