## Abstract

As the sentinels of the immune system, dendritic cells (DC) play a critical role in initiating and maintaining appropriate T cell responses through capture and presentation of antigen, costimulation, and mediator release. Although much is known about certain aspects of DC function, the exact relationship between lung epithelial DC precursor populations in the blood and their functional role in antigen presentation are not clearly understood. I created an age-structured mathematical model for DC trafficking into the lung to address this question. While capturing experimentally observed system dynamics, I found that blood DC are preferentially recruited over blood monocytes. For short-lived antigens, the model results suggest that lung epithelial DC derived from blood DC exhibit a 625% increase in antigen density compared with those derived from blood monocytes. Finally, these results motivate future experimental studies to clarify aspects of DC trafficking in the lung.

- mathematical model
- antigen presentation
- immunosurveillance

as the sentinels of the immune system, dendritic cells (DC) play an important role in initiating and maintaining T cell responses through antigen capture and presentation, costimulation, and mediator release (2, 25, 27). Human biopsy data suggest that the majority of DC in the lung epithelium are myeloid in origin, derived from either blood monocytes (BM) or blood dendritic cells (BD) (24). Although these DC precursor cells can be easily assayed in the blood, their relative contributions to the DC population within the lung epithelium and their functional roles in driving an immune response are unknown. Theoretical models are invaluable for exploring possible hypotheses, given the inherent difficulty in studying the human immune system (66).

Most theoretical studies of the immune system have focused on T cell differentiation with a minimal representation of antigen presentation (10, 11, 32, 39). Despite the key role played by DC in regulating immunity, the theoretical framework for DC dynamics has received little attention.

In previous mathematical treatments of DC, DC were simply vehicles for transporting antigen to the lymph nodes. After arrival to the lymph nodes, DC present antigen peptides to T cells that express corresponding T cell receptors. The immune response is then controlled by selecting an appropriate repertoire of non-self-reactive T cells and eliminating self-reactive T cells. As antigen transport vehicles, the total number of DC is typically modeled. In addition, previous mathematical models for DC in the lung treat the DC population as a continuum whereby the behavior of a particular cell is described by the population average (32).

Alternatively, the context in which DC are exposed to antigen plays an important role in shaping a resulting immune response. Therefore, a mathematical model for DC must not only capture the ability to transport antigen but also provide a way to convey the context in the peripheral tissue under which the antigen was captured. It has also been shown that the sensitivity of DC to environmental stimuli is highly dependent on the maturational age (4, 7, 63).

In contrast to continuum models, physiologically structured models provide a framework for representing the interplay between subpopulations with unique individual characteristics. Dynamically tracking these functionally unique subpopulations records the environmental context in which antigen is captured. Thus, given the dynamic nature of the DC population, the appropriate computational paradigm for representing DC populations is a model structured by maturational age.

In the following sections, I will develop a model for DC trafficking in the human lung epithelium using an age-structured framework. The resulting model will be calibrated to and validated against existing experimental data. Finally, I will illustrate how an age-structured model can be used to partition cellular behavior based on maturational age. Specifically, I will partition antigen uptake by DC subpopulations to predict the functional roles of the different DC precursors on antigen presentation in the lymph nodes.

## METHODS

### Model Formulation

This modeling effort focuses on representing the dynamics of the myeloid-derived DC population in the lung epithelium, as myeloid DC play a major role in the human lung under chronic and acute inflammatory conditions. A model for DC population maintenance and enhancement in response to inflammatory stimuli in the lung epithelium is composed of the regulated transport of three cell populations [BD, BM, and lung dendritic cells (LD)] and between three compartments (the blood, lung epithelium, and lymphatic system). In the following sections, I will develop a mathematical model that captures the dynamics of these cell populations.

BD and BM are recruited into the lung epithelium compartment to become LD. In general, the following population balance equation represents the dynamics of the two blood populations: (1) The “other losses” term represents the loss in the blood cell population due to recruitment to other tissues or cell death. With the use of a pseudo-steady-state approximation, the other losses term is replaced by the steady-state synthesis and lung recruitment rates. The resulting ordinary differential equation (ODE) for the BM population is then (2) where BM_{syn} is the synthesis rate of new BM, RS(*t*) is the instantaneous recruitment signal, *k*_{BM} is the sensitivity of BM to the levels of the recruitment signal, and RS(*t*_{SS}) is the steady-state lung recruitment signal. Similarly, for the BD population, the ODE is given by (3) Although recruited monocytes can potentially mature into macrophages depending on the tissue microenvironment (28, 36, 40, 41), I assume that the sole fate of recruited monocytes is to become LD. When the equations for BM and BD are combined with the equations for LD, the scaling factors LD_{tot}, BM_{tot}, and BD_{tot} are introduced to account for the different compartment sizes. These scaling factors are defined as the total number of LD, BM, and BD, respectively.

On the basis of the rapid recruitment of the blood populations into the tissue, presynthesized proteins, like P selectin, play important roles in cellular recruitment in response to inflammatory stimuli as they have expression kinetics capable of initiating the recruitment process. Presynthesized P selectin is stored in Weibel-Palade bodies in endothelial cells and can be rapidly mobilized to the cell surface within minutes on stimulation. In a cultured human umbilical vein endothelial cell system, P-selectin expression was observed to increase 2.3-fold after histamine stimulation (12, 17, 29, 51). After ∼5–10 min, surface expression reaches a maximum. Disappearance of the surface P selectin via endocytosis has been observed to occur after 30–60 min. The equation for RS(*t*) captures these acute kinetic features and is described by (4) where *k*_{basal} and *k*_{stim} are the basal and stimulated accumulation rates for the recruitment signal. The two error functions (erf) approximate a square wave starting at *t*_{RSon} and ending at *t*_{RSoff}. The rate of disappearance of the recruitment signal is given by the product of RS(*t*) with the signal decay rate constant *k*_{decay}. The steady-state value for the recruitment signal is equal to one.

Whereas BM and BD are represented as an averaged population, the LD population requires a more complex representation to capture its dynamic characteristics. The equations for LD must account for the unique characteristics of the individuals in the population. Physiologically structured models integrate mechanistic information at the individual level to give predictions of the behavior at the population level. The characteristics of each individual in a population are described by a set of physiological variables representing the internal and external characteristics of the individual. The characteristics of DC include CD antigen expression, state of maturation, and antigen loading. The state of an individual can be defined as the particular values of the physiological variables.

A relatively simple two-variable physiologically structured model is used here, with time and maturational age being the two independent variables. Additional complexity can be included, resulting in the formulation of multivariable models and, correspondingly, a more complex mathematical treatment is required. Yet, the two-variable model introduces another independent variable that describes the cell population besides time. In the case of DC, the maturational age is the second independent variable that is used to differentiate between cells of the population.

CD antigen expression has been historically used to characterize the maturational age of a DC and identify it with a particular subset (48). These external characteristics vary dynamically, as summarized in Fig. 1. The maturational age is obtained by a particular combination in expression of these CD antigens, assuming that the CD antigens do not vary independently. Thus the set of external characteristics can be replaced by a variable representing maturational age (*a*). The variable *a* ranges between 0 and 1, where 0 corresponds to a monocyte and 1 corresponds to a fully mature DC.

Using an age-structured mathematical framework, I have developed a specific model for DC trafficking in the epithelial lung tissue composed of eight equal maturational age-based subpopulations. A schematic diagram of this model is shown in Fig. 2. Derivation of the detailed equations is described in the appendix. A summary of the model variables is found in Table 1. Parameter values for each of the eight equal subpopulations are found in Table 2.

With the use of this age-structured model of the DC population, the overall behavior of lung DC can be partitioned based on the functional ability of the different subpopulations. Partitioning of antigen uptake onto the DC subsets illustrates the unique capability of an age-structured model in connecting individual- to population-based behavior. As highlighted in Fig. 3, the ability of DC to process antigen is also highly dependent on maturational age. A schematic diagram of the model for dynamically tracking internalized antigen is shown in Fig. 4. A contextual record of the lung epithelial milieu is represented by the inheritance of internalized antigen from the parent subset (DC_{i−1}) and to the child subset (DC_{i+1}). The rate of inheritance is matched to the rate of DC maturation by *k*_{mat}. Each DC subset can also internalize antigen based on its propensity and the prevailing environmental context, such as the antigen concentration [AG(*t*)]. Anantigen uptake rate [*k*_{uptake}(*ā*_{i})] can be assigned to each DC_{i} based on FITC-labeled dextran uptake reported as a function of maturational age (4). Values for *k*_{uptake}(*ā*_{i}) for each DC_{i} are shown in Table 3. The derivation of the equations for antigen uptake and tracking for each of the DC_{i} is also described in the appendix.

Parameter values were chosen based on appropriate measurements in humans, if available. However, in the absence of specific kinetic measurements in humans, data obtained from rat studies were used. The estimates for unknown parameter values (BD_{syn}, BM_{syn}, *k*_{mat}, *k*_{apop}, *k*_{BD}, *k*_{BM}, and LD_{tot}) were determined using the set ODEs and appropriate experimental data, guided by the results obtained after parameter identification. Summed squared error between experimental and simulated measurements was used to determine goodness of fit.

### Parameter Identification

An important aspect of sensitivity analysis is to identify model parameters that can be uniquely determined from the available data. A priori identifiability is one approach to select identifiable parameters given unlimited information about the modeled system (1, 23). This assessment of parameter identifiability is intimately linked to the sensitivity of the model results to changes in parameter values. Local sensitivity analysis is defined as the relative change in model variables (*y*_{j}) in response to a perturbation of the model parameters (*x*_{i}) around a local optimum (45). The resulting sensitivity function is defined as (5) where the partial derivative values are scaled by the parameter value and the maximum value of the model variable during the simulation. In practice, the sensitivity function is approximated by obtaining the sensitivity measure S_{ij} at a set of discrete time points, *t*_{k}. In addition, measurements of only a subset of the model variables may be obtained experimentally. A reduced sensitivity matrix (M) is then constructed using a subset of the sensitivity measures [*S*_{iĵ}(*t*_{k})] such that (6) A set of correlation coefficients between model parameters is calculated from M. Parameters that are locally identifiable have correlations with all other parameters between −1 and 1. Parameters that are not locally identifiable, termed a priori unidentifiable, have correlations of exactly +1 or −1 with at least one other parameter.

## RESULTS

### Simulating the Blood Population Dynamics

Before the response to an acute stimulus was simulated, parameter values and rate constants must be determined for steady-state conditions. Some of the most readily available data are measurements of the blood populations. BM comprise ∼5% of the cells in human blood with a median density of 5.0 × 10^{8} cells/l (42). In contrast, BD are much fewer in number than BM. From blood samples of patients with asthma, Upham et al. reported the density of blood DC (HLA-DR^{+}, CD33^{+}, CD14^{−}, and CD16^{−} cells) to be 17.4 ± 5.4 *×* 10^{6} cells/l. Other groups have estimated a normal range for BD to be from 0.15% to 0.70% mononuclear cells or 3 × 10^{6} to 17 × 10^{6} cells/l (8, 9, 19). With the assumption that the average human has 5.6 liters of blood, the value for the total population of BM and BD is 2.80 × 10^{9} cells and 9.74 × 10^{7} cells, respectively.

The half-lives of BM have been studied by numerous groups (16, 35, 56–58, 62). However, many of these studies were conducted with a single dose of radiolabeled tracer. Van Furth (56) stated that cell transit times in humans were similar to those observed in rats and much slower than shown in mice. Thus the continuous radiolabeling study in rats reported by Whitelaw (61) provides an estimate of the human BD and BM turnover rate. Under normal steady-state conditions, the syntheses of the blood populations are exactly balanced by the rates of cellular migration into the tissue compartments. Thus the blood half-life and synthesis rate of BM were determined to be equal to 77.9 h and 2.49 × 10^{7} cells/h, respectively. When the synthesis rate for BM is set to zero, the mathematical model is compared against the experimental data of Whitelaw in Fig. 5. A similar estimate of the half-life of BM has been reported in humans using a pulse-labeling technique (62). Because both the BM and BD come from a common precursor population, the cellular BD turnover rate was assumed to be the same as for BM, resulting in a BD synthesis rate equal to 9.13 *×* 10^{5} cells/h. A summary of parameter values is shown in Table 2.

### Simulating the LD Dynamics

#### Steady-state turnover.

Various experimental groups have tried to measure the dynamic responses of the LD populations to various experimental conditions (21, 26, 50). Much of the understanding of the steady-state turnover and rapid influx after acute challenge is based on a series of studies from the Holt group (20, 24, 33, 34). To determine the steady-state turnover of DC populations in the rat lung, Holt et al. determined that the resident population of DC in the epithelium declined by 85% in 72 h after x-irradiation of bone marrow renewal populations (20). Although these studies were performed in a rat system, similar kinetics are expected in humans (24, 62). With these data, values for *k*_{mat} and *k*_{apop} were determined to be 0.209 h^{−1} and 0.163 h^{−1}, respectively. In an alternative approach to the elimination of bone marrow renewal populations, they also found a similar reduction in DC populations after dexamethasone administration. As shown in Fig. 5, simulating the experimental protocol described by Holt et al. (20) gives good agreement with the experimental findings.

#### Acute challenge response.

In characterizing an acute challenge response, McWilliam et al. (34) measured the kinetics of DC migration into rat tracheal epithelium after 60-min inhalation exposure to heat-killed bacteria. Relative to background, they observed a 3.5-fold increase peaking at 1 h after cessation of antigen exposure. Similarly, Gong et al. (15) also found a two- to threefold rise in airway DC after IFN-γ instillation. Similar behaviors were observed in ovalbumin-sensitized (33) and Hen-egg lysozyme-immunized rats (64). Similarly, human biopsy data obtained 4 h after an allergen challenge revealed a 3.5-fold increase in LD density relative to prechallenge levels (24).

Upham et al. (54) recently reported the changes in blood DC precursors after allergen challenge. These kinetic measurements provide additional insight into the dynamics of DC recruitment. At 3 h after allergen challenge, Upham et al. observed a 37% drop in BD and a 5% drop in BM. After 24 h, the BD concentration had returned to 78% of the prechallenge level, whereas the BM concentration had returned to normal.

The final step is the estimation of the remaining parameter values for LD_{tot} and the recruitment rate constants (or sensitivities) for BM and BD in response to RS(*t*). The recruitment rate constants are represented by the parameters *k*_{BM} and *k*_{BD}. The values for these parameters obtained using the kinetic data for an acute challenge protocol are shown in Table 2.

Simulation results for both the BD and LD populations in response to the acute challenge protocol are in good agreement with the in vivo studies reported by the Holt group (34, 54), as shown in Fig. 6. The simulation results also suggest that there is an enrichment of BD-derived DC in the lung compared to BM-derived DC. As shown in Fig. 7, the simulation suggests that ∼20% of the LD originates from BD.

#### Parameter identification.

Given the different protocols used in the calibration data, the antigen challenge protocol was used to identify the parameters because it provides the best estimate of the parameter sensitivities. The sensitivity matrix was also reduced to include only sensitivities for the variables BM, BD, and LD because additional variables are not presently reported. In addition, a subset of parameters was analyzed to highlight the model parameters determined in this study. The correlation coefficients for the parameters are shown in Table 4, with significant values (values > 0.965 and *P* < 0.05) highlighted in bold. The results suggest that *k*_{BM} and LD_{tot} are correlated and thus not identifiable. However, the maximum sensitivities, shown in Table 5, suggest that the different experimental protocols (irradiation vs. acute challenge) used in calibration enable these parameters to be uniquely identified. Specifically, *k*_{BM} can be determined based on the dynamic response of BM to the lethal irradiation and antigen challenge experiments, as the maximum sensitivity of LD_{tot} on BM and BD is equal to 0.

#### Antigen processing roles for BM and BD.

One of the motivations for this study is whether there is a functional consequence for differential recruitment of BM vs. BD. To assess this question, the antigen uptake capacity was partitioned into the eight DC subsets of the trafficking model. The antigen densities in the mature DC subset (DC_{8}) derived from either a BM or BD are measured separately, using equations that dynamically track internalized antigen in the DC population. The acute challenge protocol was used in the simulation. Because antigens display a range of half-lives in the lung, the antigen density is reported over a range of antigen decay rate constants (*k*_{Agdeg}). For BM-derived mature DC, Fig. 8, *top*, intuitively shows that, as the antigen decay rate decreases, the antigen density increases. In Fig. 8, *bottom*, I report the increase in antigen density in BD-derived DC relative to BM-derived DC. At 50 h after the acute challenge, BD-derived DC exhibit a 625% increase in antigen density compared with BM-derived DC for an antigen decay rate constant of 0.6 h^{−1}.

## DISCUSSION

Although physiologically structured models have been proposed since the mid-1960s (13), they are seldom used because of the difficulty in obtaining appropriate experimental data and the models' mathematical complexities. In particular, physiologically structured models require information at the single-cell level, such as single-cell state-to-state transition rates and partitioning functions. Physiologically structured models have been presented for blood cell populations (3, 31), but this is the first example of applying an age-structured framework to a tissue population. The initial steps in developing such a model, and the objectives of this study, are to integrate available knowledge about DC trafficking into a computational framework, describe how state-to-state transition rates can be obtained experimentally, and provide testable hypotheses. In the future, recent advances in flow cytometry (18, 22) and multicolor immunohistology (53) may provide the requisite data to further validate an age-structured model of DC trafficking in the lung.

Qualitatively, McWilliam et al. (34) observed an age structure in the lung epithelial DC population, stating “at the 8-h time point, the morphology of the intraepithelial Ia^{+} cell population had changed markedly towards a more pleomorphic form reminiscent of DC in their 'veiled cell' form and, by 24 h, many of these cells appeared to have developed into typical mature, highly branched DC.” By simulating the LD population as age structured, we see that DC are equally distributed among all maturational states under steady-state conditions. However, after antigen exposure, a large bolus of immature DC enters the lung epithelium, matures in unison, and exits the tissue by emigrating into the lymphatic system. This bolus of DC dominates the population characteristics and gives rise to the observation by McWilliam et al. Our results suggest that the average transit time of a BM-derived DC through the lung epithelium is ∼24 h, significantly quicker than may be inferred from in vitro studies.

One of the objectives of this modeling study is to gain insight into the biological system by understanding the source of systematic deviations between simulated and experimental data. Testable hypotheses are obtained by considering these systematic deviations in the context of the model structure. In the following sections, we examine in more detail the two different trafficking scenarios (steady-state and acute challenge protocols) used in developing this model.

In general, there is good agreement between the simulated and experimental protocols for assessing the steady-state turnover rates of the blood precursor populations and the LD population. Thus the hypothesis that BM and BD exhibit similar turnover rates and contribute to the lung DC pool is consistent with the experimental data. In examining the experimental data more closely, the disappearance of LD reported by Holt et al. (20) is at a slightly higher rate than the rate of BM disappearance measured by Whitelaw. Because it is hypothesized that LD are derived from blood precursors rather than self-replication, it is unexpected that LD decline at a higher rate than their precursor populations. One possible explanation could be that the labeling of cells by Whitelaw et al. may have been inefficient, resulting in a slower decline in unlabeled BM. This suggests that the true half-life of BM may be lower than can be estimated from Whitelaw's data and may actually be more consistent with the observations of Holt et al. (20). An alternative explanation is that the lethal irradiation protocol may have increased the *k*_{mat}, resulting in a lower apparent half-life than would be expected under steady-state conditions. Given the uncertainty in the turnover rates of BM and BD, it may be insightful to apply more modern immunohistochemical and radiological techniques to measure the half-lives of these individual populations in the blood and tissues.

For the acute challenge protocol, there is good initial agreement between the simulation and experimental results. At the 6- and 24-h time points, the simulation suggests that the blood populations recover to their prechallenge levels through additional secondary feedback pathways. Hematopoietic regulation by inflammatory cytokines that become elevated in the systemic circulation could be a possible feedback mechanism. In addition, the recruitment signal is comprised solely of an early event. However, RS may also include a delayed contribution due to the time lags introduced by protein synthesis (signal transduction, DNA transcription, RNA translation, and protein release) in response to the antigen challenge. Inclusion of a delayed RS and improved regulation of BM and BD would both increase the LD population. Yet animal and human data suggest that the LD population stabilizes at an elevated level shortly after the antigen challenge. McWilliam et al. (34) also observed that the LD density remained at an elevated level for ∼24 h before slowly declining to baseline at 72 h. Consequently, this suggests that *k*_{mat} is not constant but changes dynamically.

On a percent basis, I estimated that BD comprise only 3.4% of the total blood precursors that can become LD. In contrast, the results from this mathematical model for DC trafficking suggests that 20% of the LD that accumulates in response to an inflammatory stimulus originates from BD. This selective enhancement of BD in the tissue implies that, on a per-cell basis, BD are preferentially recruited into the tissue over monocytes. The preferential recruitment of BD is quantified in the model by a nine times larger value for *k*_{BD} than for *k*_{BM}, where *k*_{BD} and *k*_{BM} are the sensitivities to the recruitment stimulus of BD and BM, respectively. One possible explanation could be the upregulated expression of chemokine receptors or adhesion molecules in BD vs. BM. Vanbervliet et al. (59) recently reported the differential expression of chemokine receptors on DC precursor populations. Thus the calculated enrichment of BD in the tissue may be a result of the upregulation of CCR2, CCR5, CCR6, or CXCR4 on BDC compared with BM.

As shown in Fig. 8, the consequences of BD vs. BM recruitment may be profound. As described earlier, the functional relationship between maturational state and the ability of DC to process and present antigen is based on the experimental studies of various groups (4, 14, 44, 65). Because BD could be considered more functionally mature than BM, it follows that BD could have a more developed ability to process antigen. Thus newly recruited BD may provide rapid sampling of short-lived antigens in the lung. In contrast, BM-derived DC are unable to sample adequately short-lived antigens as they are cleared from the tissue before BM-derived DC acquire the full ability to process antigen. Because long-lived antigens may be difficult to distinguish from self-antigens, BM-derived LD may play a tolerogenic role in controlling immunity in the lung.

Another important value determined from modeling DC trafficking is the size of the LD pool (LD_{tot}). The value obtained in this study can be cross-validated against alternative experimental approaches. With the use of immunohistochemical approaches, various groups have estimated the density of the epithelial DC network and also found that the density of DC is strongly dependent on the location in the tissue and the airways (20, 33, 34, 46). Holt and coworkers (20, 33, 34) reported that, in a rat system, the DC epithelial network is composed of between 500 and 800 DC/mm^{2} of epithelial surface. In an earlier study, Schon-Hegrad et al. (46) reported that the density of DC varies from 600–800 per mm^{2} epithelial surface in the large airways to 75 per mm^{2} in the epithelium of the small airways of the peripheral lung. With these values, the density of LD (in cells/mm^{2}) as a function of generation (*Z*) is defined using a Weibel representation of the airway structure (60): (7) The total number of DC in the conducting airways of a adult human lung can be estimated to be equal to 1.8 × 10^{8} DC/lung by combining *Eq. 7* with an estimate of the epithelial surface area of the airways as a function of generation (60).

Using an alternative approach to measure the size of the LD population, Gong et al. (15) used microdissection techniques to isolate DC from a rat lung epithelium and phenotypically compared them to DC isolated from the rat parenchyma. In the digests from the rat lung, they found 6.2 × 10^{3} epithelial DC/lung and 6.8 × 10^{4} parenchymal DC/lung. With lung digestion techniques, a 50% yield of LD is typically expected. These densities can be scaled to humans by assuming that lung size is proportional to body weight and that a typical weight for a rat and human is 0.17 and 70 kg, respectively. This microdissection approach suggests that the human LD pool is 6.0 × 10^{7} DC/lung.

As shown in Table 2, LD_{tot} is estimated to be equal to 8.66 × 10^{7} DC/lung and agrees with both of these estimates. This suggests that, in addition to a present hypothesis where large-scale migration of DC from immunologically perturbed tissues can be attributed in part to the mobilization of interstitial and epidermal DC (30, 43), BD provide a significant source for the observed increases in LD populations in response to inflammatory stimuli.

In conclusion, the primary focus of this model for DC trafficking in the lung epithelium is to understand the contributions of myeloid blood precursor populations to homeostasis and inflammatory response in the human lung. I estimated the contributions and possible mechanistic roles for BM and BD under conditions of acute and chronic inflammation by integrating experimental observations into a quantitative age-structured framework.

The new age-structured mathematical model presented here for the trafficking of LD captures experimentally observed system dynamics and provides insight into some of the governing biological phenomena. In particular, presynthesized proteins that facilitate cellular recruitment, like P selectin, may play a significant role in the recruitment of DC in response to acute inflammatory stimulus based on the kinetics of expression. These results also suggest that LD are primarily recruited from blood precursor populations, where BD are preferentially recruited over BM. Because BD display enhanced antigen processing capabilities, BD-derived mature LD provide a 625% increase in antigen density for short-lived antigen compared with BM-derived DC. As a result, BD-derived LD may play a more significant functional role in activating an adaptive immune response.

Although experimental validation of the age structure to the LD population is presently qualitative, dynamic immunohistochemical characterization based on cell surface marker expression may provide insight into the relative contribution of the blood precursor populations to LD. Quantitative characterization of the LD pool would also validate the entry point of BD into the age-structured population. Finally, dynamic measurement of the cell surface marker expression would characterize the size of each DC_{i} pool and the state-to-state transition rates.

As highlighted by this study, better understanding is needed of the differences between BM and BD in terms of recruitment, antigen processing, and hematopoietic regulation. Quantifying both BM and BD turnover would aid in understanding the hematopoietic regulation of these populations. In addition, the feedback on hematopoietic regulation of blood populations requires further elucidation. Because most of the data used in the calibration of the turnover rates were derived from older experimental techniques, it may be insightful to apply more modern immunohistochemical and radiological techniques to the study of the half-lives of these DC populations in the blood and tissues.

Ultimately, understanding the dynamics of DC trafficking may lead to therapeutic strategies to skew the immune response in a desired direction (38). Yet, current approaches for developing DC-mediated anti-cancer therapies have failed to bear fruit (5, 47, 52). Given the complex multivariate relationship between environmental inputs and DC-controlled outputs, creation of a cell population-based model of this process to explore alternative multivariable control approaches (37) for therapeutic advantage would be invaluable. This study is an essential first step toward that goal.

## APPENDIX

### Derivation of Model Equations

The following section describes the derivation of the equations used in this model. The total DC population in the lung can be described as a continuous distribution of cells, assuming that there are a statistically significant number of DC in the lung epithelium. The population of DC in the lung epithelium at a given point in time is then defined in terms of the density of DC with maturational ages between *a* and *a*+*da* [DC(*t*,*a*)*da*]. Furthermore, the change in DC density is expressed as a function of the recruitment from the blood *R*(*t*,*a*) and emigration rate into the lymphatics [*E*(*a*)] such that (A1) Because the density of DC depends on both time and maturational age, the derivative of the DC density is equal to (A2) By combining*Eqs. A1* and *A2* and dividing by *dt*, the change in DC density is represented by the partial differential equation: (A3) *k*_{mat} replaces the fraction *da/dt* by assuming that the maturation rate is constant.

The total density of LD is obtained by integrating over maturational age (A4) The partial differential equation in *Eq. A3* is converted into a collection of coupled ODEs by discretizing LD into *n* subpopulations such that (A5) where each DC_{i} is defined by integrating *Eq. A4* over a range of ages such that (A6) The maturational age of each DC_{i} is represented by the average age (*ā _{i}*) equal to (

*a*

_{i−1}+

*a*

_{i}). Values of

*ā*

_{i}are shown in Table 3. The generic ODE for LD subpopulation

*i*is described by (A7) A set of similar coupled ODEs represents the entire lung epithelial DC population. However, the first and

*n*th subpopulations require special treatment. The equation for DC

_{1}reduces to (A8) as the “flow” into this subpopulation occurs solely through the R(

*t*,

*a*). Similarly, the equation for DC

_{n}simplifies to (A9) as the “flow” out of this subpopulation does not occur based on maturation but programmed cell death. The rate of programmed cell death (apoptosis) is represented by the parameter

*k*

_{apop}.

DC may produce mediators that contribute to the inflammatory milieu, resulting in nonlinear feedback through the DC recruitment term R(*t*,*a*). However, because DC are not the only contributors to the inflammatory response, I assume that R(*t*,*a*) is independent of DC(*t*,*a*). Thus the recruitment term is expressed as (A10) where BP*ā*_{i}(*t*) is the blood precursor population that gets recruited into DC_{i}, and Sens(*ā _{i}*) is the sensitivity of BPā

_{i}(

*t*) to the recruitment stimuli. Values for Sens(

*ā*) are shown in Table 3.

_{i}As discussed in the text, RS(*t*) represents the combined dynamic effect of the introduction and elimination of an environmental insult, release and elimination of presynthesized mediators like histamine from mast cells, and mobilization and endocytosis of P selectin on the endothelial surface. Although the recruitment process is generally well characterized, it is difficult to estimate the magnitude of this recruitment signal. For simplicity, the steady-state value of RS(*t*) was assigned a value of 1. The value for *k*_{decay} was chosen to have a half-life of ∼15 min, similar to that observed for P selectin. Once *k*_{decay} was assigned, *k*_{basal} can be obtained by assuming a steady-state value of RS(*t*) of one. The two-error function terms represent a square-wave source of magnitude equal to one starting at *t*_{RSon} and ending at *t*_{RSoff}. Numerically, the G_{RS} terms “soften” the onset and cessation of the source term such that the square wave is a continuously differentiable function. The resulting profile for RS(*t*) is shown in Fig. A9.

The emigration rate simplifies to (A11) where *E*(*ā _{i}*) is the average emigration probability for DC

_{i}. The average emigration probability is derived from the degree of downregulation of tissue homing and upregulation of lymphatic homing chemokine receptors as a function of maturational age (6, 49). The values for

*E*(

*ā*)are shown in Table 3. As summarized in Fig. 2, the particular model equations used in the simulation for BM, BD, and DC

_{i}_{i}are consistent with this general formulation of an age-structured framework. The equations include specific adaptations to account for DC trafficking to the lymph and are as follows: (A12) (A13) (A14) (A15) (A16) (A17) (A18) (A19) (A20) (A21) For the acute challenge protocol, the concentration of antigen is described by (A22) where the error functions approximate a step function of height 1 starting at

*t*

_{Agon}and ending at

*t*

_{Agoff}. The rate constant

*k*

_{Agdeg}parameterizes the decay rate in antigen concentration. Values for the parameters are shown in Table 2.

Following from Fig. 4, the schematic diagram for tracking a cumulative characteristic of DC_{i}, antigen uptake and tracking for each of the DC_{i} are given by the following equations: (A23) (A24) (A25) (A26) (A27) (A28) (A29) (A30)

### Parameter Identification

As mentioned in the text, an important aspect of sensitivity analysis is to identify model parameters that can be uniquely determined from the available data. This concept is commonly referred to as the identifiability problem and precedes the estimation of parameter values, given a model of a system and specific input-output measurements. A priori identifiability is one approach to select identifiable parameters given unlimited information about the modeled system (1, 23). This assessment of parameter identifiability is derived from the sensitivity of the model results to incremental changes in parameter values. The sensitivity of the model results are quantified by a local sensitivity analysis. A local sensitivity analysis is defined as the relative change in model variables (*y _{j}*) in response to a “local” perturbation of the model parameters (

*x*) around a particular set of parameter values. In this case, the set of parameter values used were the values that exhibited the best fit between the model results and experimental data. Although identifiability analysis can be performed at any point in parameter space, it is typically measured at an optimum (45). Although values of the sensitivity measure, S

_{i}_{ij}, are determined at a particular point in time, a sensitivity function quantifies the time-dependent changes in a particular dependent variable given perturbations in a particular parameter value. The resulting sensitivity function is defined as (A31) where the partial derivative values are scaled by the parameter value and the maximum value of the model variable during the simulation. This scaling normalizes the sensitivity function such that the results are not biased by the relative values of the variables and parameters due to an arbitrary selection of units. In practice, the sensitivity function is approximated by obtaining S

_{ij}at a set of discrete time points,

*t*

_{k}. Furthermore, measurements of only a subset of the model variables may be obtained experimentally. Thus a reduced sensitivity matrix (M) can be constructed using a subset of the sensitivity measures [S

_{iĵ}(

*t*

_{k})] such that (A32) In this sensitivity matrix, the columns correspond to the sensitivity measure for the different variables

*y*

_{j}for a particular parameter

*x*

_{i}. Generally, correlation statistics can be calculated between independent variables (parameters) given a matrix whose rows are the observations and columns are the independent parameters. In the simulation tool Matlab, this is done with the command “corrcoef.” A set of correlation coefficients between model parameters is calculated from M. Parameters that are locally identifiable have correlations with all other parameters between −1 and 1. A correlation of −1 means that the parameters are inversely related to each other in their effects on the variables. Conversely, a correlation of 1 means that the parameters are directly related to each other in their effects on the variables. Parameters that are not locally identifiable, termed a priori unidentifiable, have correlations of exactly +1 or −1 with at least one other parameter. In other words, if the perturbation of two parameters result in the same set of observations (and thus are highly correlated), it will be difficult to determine the parameter values uniquely from a set of experimental observations. To address situations where the model parameters are unidentifiable a priori, one can simplify the model to eliminate redundant parameters or obtain data that discriminates between the parameters. In calibrating the model, experimental data were used that discriminate between the parameters that might be considered unidentifiable for the antigen challenge simulation.

## Footnotes

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 © 2006 the American Physiological Society