M Tiphine

B Vezzoni

F Alvarez-Velarde

B Carlier

D Wojtaszek

G Glinatsis

A Brolly

B Dixon

R Eschbach

C Coquelet-Pascal

D Freynet

G Edwards

B Hyland

S Häkkinen

T Viitanen

F Gabrielli

K Ono

S Cornet

The Effects of the Uncertainty of Input Parameters on Nuclear Fuel Cycle Scenario Studies

Nuclear systems, composed of reactors and fuel cycle facilities, are complex and in constant evolution. Fuel cycle scenario studies aim at evaluating the consequences of decisions taken in the near term on the eventual outcomes of a chosen fuel cycle.

Scenario studies require a large amount of data and hypotheses, such as the energy demand and its evolution or the performances and characteristics of existing and future technologies. Each hypothesis, represented by a parameter, is a potential source of uncertainty that can impact the scenario outcomes. The Nuclear Energy Agency's Expert Group on Advanced Fuel Cycle Scenarios (NEA/AFCS) recently completed an activity that used sensitivity studies to evaluate the impact of uncertainty of the primary scenario input parameters on scenario results.

At the beginning of the study, a base case scenario was specified and evaluated. It considers a PWR UOX reactor fleet using low enriched uranium oxide fuel and its gradual replacement by a sodium fast reactor (SFR) fleet using uranium/plutonium mixed oxide fuel. Each member of the group evaluated this scenario with its own scenario code (COSI, COSAC, FAMILY, SITON, TR_EVOL or VISION). It is noteworthy that the different codes show a good consistency in the results and that the differences have been investigated and documented in order to ensure a strong basis for comparison in the sensitivity study.

Then, 17 key input parameters and their range of variation as well as 22 scenario outputs were identified and sensitivity analyses were conducted. The key parameters were chosen amongst the general assumptions on the scenario (energy demand, introduction date of fast reactor…), reactors and facilities characteristics (burn-ups, losses at reprocessing…) and some minor actinides recycling parameters. As far as possible, each parameter was varied independently from the others. This sensitivity study shows the relative impact of each parameter on each output metric.

In order to compare the large amount of results, two methods of summarization were adopted, as described in the paper. These comparisons have shown that some parameters, such as the energetic production or the introduction date of fast reactors, have a strong impact on the results whereas others, like the UOX fabrication time, have almost no impact on the studied scenario.

Introduction

Under the auspices of the NEA Nuclear Science Committee (NSC), the Working Party on the Scientific Issues of the Fuel Cycle (WPFC) has been established to coordinate scientific activities regarding various existing and advanced nuclear fuel cycles, including advanced reactor systems, associated chemistry and flowsheets, development and performance of fuels and materials, accelerators and spallation targets. Various Expert Groups were established to cover the abovementioned topics.

The Expert Group on Advanced Fuel Cycle Scenarios (EGAFCS) was created in 2010, replacing the Expert Group on Fuel Cycle Transition Scenarios Studies, to study R&D needs associated with the transition from current or future advanced nuclear fuel cycles. The objectives of the expert Group are to assemble, organise and understand the scientific issues of advanced fuel cycles, and to provide a framework for assessing specific national needs related to the implementation of advanced fuel cycles.

After conducting a benchmark study to compare existing codes in terms of capabilities, modelling and results [START_REF]Benchmark Study on Nuclear Fuel Cycle Transition Scenarios Analysis Codes[END_REF], the expert group performed a benchmark study to identify and communicate the impact of uncertainties in fuel cycle analyses using systems codes [2]. The purpose of this study was to assess and quantify the importance of input parameters in scenario analysis. The benchmark was conducted in different phases: definition of a base case scenario, system codes normalisation, parameters identification, parametric studies, and uncertainties representation.

Description of the study

Fuel cycle scenario analysis is a common method for identifying and communicating potential nuclear energy futures, especially when assessing the impact of new technologies. System codes are used to assess the behaviour of different nuclear energy systems and to examine scenarios involving transition of infrastructure (reactors, fuel cycle facilities, etc.) to accommodate growth and adopt new fuel cycle management approaches (e.g. recycling).

By their nature, scenario analyses involve assumptions about the future, and the systems codes used in scenario analyses are imperfect representations of the real world. Uncertainties inherent in scenario specifications and system codes impact the accuracy of analyses. There are several sources of uncertainty in fuel cycle analyses using system codes, including scenario assumptions, facility operational values, reactor core physics calculations, etc. and their impacts vary based on the fuel cycle and the performance metrics of interest. The purpose of this study is to systematically identify these sources of uncertainty and use sensitivity studies to assess their impacts on system level results.

Scenario codes

The study involved the use of multiple system codes, some of them used by more than one organisation. The different codes used to obtain those results are the following:

-COSI, used by CEA, ENEA and KIT [START_REF] Coquelet-Pascal | COSI6: a tool for nuclear transition scenarios studies and application to SFR deployment scenarios with minor actinide transmutation[END_REF]; -COSAC, used by AREVA [START_REF] Carlier | COSAC: a computational code dedicated to fuel cycle back-end simulations in AREVA[END_REF];

-FAMILY 21, used by JAEA; -SITON, used by EK [START_REF] Brolly | Development and application of SITON, a new fuel cycle simulation code[END_REF]; -TR_EVOL, used by CIEMAT [START_REF] Álvarez-Velarde | TR_EVOL, upgrading of EVOLCODE2 for transition scenario studies[END_REF]; -VISION, used by CNL and INL [START_REF] Jacobson | Verifiable Fuel Cycle Simulation Model (VISION): A Tool for Analysing Nuclear Fuel Cycle Futures[END_REF].

All of the participants modelled the base case scenario and adjusted their models to produce similar behaviour. The analyses of the base case scenario by multiple organisations uncovered any remaining ambiguities in the scenario specification and helped clarifying the scenario documentation in order to ensure consistent interpretation. The remaining discrepancies revealed the impact of differences in code architectures. When these differences occurred, their drivers were identified and discussed. Finally, the base case results for each code provided a basis for identification of the effects of uncertainty of key parameters.

Description of the base case scenario

This study was conducted in several phases. First, a base case scenario was identified and carefully defined. The Expert Group has learned from past efforts that a clear and concise scenario definition is necessary to reduce the potential for multiple interpretations of the specification. The resulting base case scenario is 200 years in duration, begins with a fleet of thermal reactors (PWRs) in equilibrium using low enriched uranium oxide fuel (UOX) and transitions to a fleet of sodium fast reactors (SFRs) using a uranium/plutonium mixed oxide fuel near mid-scenario. The SFR fleet reaches its equilibrium before the end of the scenario. The electricity production is maintained constant at 430 TWhe/y throughout (see Figure 1).

Figure 1: Base case scenario -installed capacity

In the base case scenario, the reactors lifetime is supposed infinite. The enrichment tail of 235 U is 0.25%. Fuel fabrication starts two years before the first load in each reactor and SFR fuel is made with a mix of depleted uranium and plutonium (Pu). The minimum spent fuel cooling time before reprocessing is of 5 years for PWR fuels and 2 years for SFR fuels. Time for reprocessing is assumed to be null. PWR UOX spent fuel reprocessing starts in year 35 and is fixed at 850tHM/y. SFR spent fuel is reprocessed when available and according to the need in Pu for fresh fuel fabrication. The Expert Group selected 22 of the most important scenario outputs as indicators for the study. Standard outputs include reactor numbers and types, mass flows of fuel and waste materials, and inventories of materials in storage and disposal.

Identification of the key parameters

The most important component of this study was the identification of the key parameters to be evaluated. The Expert Group drew from its collective expertise to develop the set of 17 parameters that most often are sources of uncertainty during the construction of scenario specifications. The second step in this process was to identify a variation range for each key primary parameter to be used as a basis for the sensitivity analysis. Again, experience was used, with each parameter varied over ranges typically used in past analyses (see Table 1).

Output from the base case scenario

The outputs from the base case scenario obtained by all participants were compared in order to identify the discrepancies between the different codes. It is noteworthy that the different codes show a good consistency in the results. As 22 outputs were studied, those for which a very good agreement between the different codes was observed are not detailed in this article. The remaining differences are detailed hereafter.

For the front-end cycle, we observed a really good agreement between the different codes for the evaluation of the natural uranium consumption, the enrichment need, the UOX fabrication need and the Pu flow for FR fuel fabrication. Slight discrepancies were observed on the estimation of the FR fuel fabrication need (see Figure 2). They are due to inconsistencies in the initial data as this estimation was made from the energy production for some codes (e.g. COSI6) or from the installed nuclear power for others (e.g. FAMILY).

Figure 2: Base case scenario -FR fissile fuel fabrication need

For the back-end cycle and the inventories, more discrepancies were observed on the estimation of outputs. One source of discrepancies is the modelization of the reprocessing. Indeed, in the scenario specifications, the UOX spent fuel reprocessing starts in year 35 and is fixed at 850tHM/y in order to produce enough plutonium for FR deployment. However, with the code COSAC, available spent fuel is reprocessed only when Pu is required and in function of the need in Pu, without the possibility to anticipate it. That creates discrepancies on the need in PWR UOX spent fuel reprocessing (Figure 3), and thus on the PWR UOX spent fuel and the interim reprocessed uranium storages (Figure 4). This issue is not encountered for the FR spent fuel reprocessing, as it follows the need in Pu for all codes. Another source of discrepancies is a slight difference in the breeding gain of the fast reactors considered by each organization. As SFR spent fuel is reprocessed according to the need in Pu, the amount and fissile content of Pu in spent SFR fuel directly impact their flow at reprocessing. Thus, differences on the breeding gain impact the need in SFR spent fuel reprocessing (Figure 5), and so the SFR spent fuel storage (Figure 6) and the Pu interim storage (Figure 6). The discrepancies observed on the total inventory of Pu in cycle and on the minor actinides inventory in waste are the consequences of the differences described previously (Figure 7). Overall, the study of the base case scenario by all the benchmark participants lead to consistent results between the different codes and provided a strong basis for the sensitivity study.

Sensitivity analysis

The primary activity in this study was to conduct sensitivity analyses on the key parameters in order to identify the effects and to quantify the potential impacts of sensibilities of the outputs to the inputs. To the extent possible, each parameter was varied independently, without change in any other part of the specification. During the sensitivity analyses, some scenarios "broke" and had insufficient fuel material for the SFR fleet. This was noted in the analysis and then an effort was made to "fix" the scenario. One option was to modify additional parameters to address an imbalance in the scenario. The other option was to add an external source of fuel material and to note the amount of extra material required to complete the scenario.

The sensitivity of the 22 outputs to the 17 identified parameters has been evaluated. With each fuel cycle metric being evaluated for each year of the scenario, we had to consider ways of comparing and displaying the results. Two methods were used and are detailed in this paper: the tornado diagrams and the sensitivity table.

To quantify the impact of the variation of an input parameter on the calculation, output indicators had to be defined. Three types of output indicators were retained:

-All the output parameters related to an annual flow (e.g. natural uranium consumption, enriched uranium needs, fuel fabrication and spent fuel reprocessing, …) were associated to output indicators defined as the cumulated value of the concerned output parameter;

-All the output parameters related to an inventory in an interim storage (e.g. spent fuel storage, reprocessed uranium storage, depleted uranium storage, …) were associated to output indicators defined as the maximum value of the concerned output parameter;

-All the output parameters related to a global inventory inside the cycle or the waste (e.g. plutonium inventory in the cycle, minor actinide inventory in the waste, …) were associated to output indicators defined as the final value of the concerned output parameter.

Tornado diagrams

The tornado diagram shows the correlation of a single output parameter to changes in the input parameters, one bar per input, where the direction of the bars on the diagram shows whether the correlation is positive or negative. The length of the bars shows the magnitude of the coupling between the parameters. There are two values of q for each input parameter: one for the case in which the given parameter is at its low value, and one for the case in which the parameter is at its high value.

The sensitivity value (see equation ( 1)) is used to quantify the change in an output indicator as a function of the change in the input parameter value, thus giving an estimation of the relative impact of each parameter on each output indicator. Note that the sensitivity cases that require changing more than one input parameter are excluded from this analysis.

) ( The results of the sensitivity study, represented with Tornado diagrams, are detailed hereafter.

) ( s ref ref s ref ref p p R R R p q - - = (

Front-end cycle

Figure 8: Sensitivity of front-end cycle results to each parameter

Overall, the beginning date of the SFR introduction has a strong impact on the front-end cycle results. Indeed, it is directly correlated to the period of operation of PWR and SFR. A later introduction of SFR results in a longer operation time of PWR and so on higher natural uranium (Unat) consumption and separative work unit (SWU) requirements. On the other hand, it reduces the SFR fleet operation time, resulting in reduced SFR fuel fabrication need and Pu flow to fabrication. For the same reasons but to a lesser extent, the duration of SFR introduction (their commissioning rate) impacts the results in the same way as their introduction date.

Obviously, the UOX enrichment tail impacts the SWU requirements. Due to the additional effort (SWUs) required to extract additional fissile material from lower assay tails, the relationships for tail enrichment are non-linear.

Due to the fabrication of final cores (resp. start-up), the PWR (resp. SFR) lifetime is negatively correlated to the PWR (resp. SFR) fuel fabrication need. This impact is slightly non-linear, as it is driven by the fraction of start-up cores over reactor life.

SFR fuel fabrication and Pu flow to fabrication are primarily impacted by the energy growth rate, which determines the total number of SFR required.

Finally, the SFR breeding ratio has a large, non-linear impact on the Pu flow through fuel fabrication. This is primarily due to SFRs with a negative breeding gain requiring an additional external Pu feed and impacting the fissile content and the quantity of the Pu.

Back-end cycle

Figure 9: Sensitivity of back-end cycle results to each parameter

Spent fuel storages are sensitive to the reprocessing parameter as a later start of reprocessing results in an accumulation of spent fuel, but a higher reprocessing rate leads to a faster depletion of the storage. PWR UOX spent fuel reprocessing parameters also impact the amount of interim Pu stored. It is not the case for SFR spent fuel as it is reprocessed according to the need in Pu.

The energy growth rate has a large impact on SFR spent fuel reprocessing as it determines the number of commissioned SFRs. This property is non-linear as with a higher growth rate there is relatively more SFR in cooling inventory that is not yet available for reprocessing.

As the SFR spent fuel flow at reprocessing is adapted according to the Pu need for fresh fuel fabrication, it depends on the amount of Pu available in the spent fuel and so is impacted by the breeding ratio. The breeding ratio also directly impacts the separated Pu inventory as a higher breeding ratio increases the available Pu without significantly changing the Pu needed.

As the SFR introduction date is directly linked to the operational time of reactors and to the amount of available PWR and SFR spent fuel, it highly impacts the PWR and SFR spent fuel flow at reprocessing. Separated Pu being used for SFR fuel fabrication, a later SFR introduction leads to an accumulation of separated Pu.

Finally, the reactors lifetimes are negatively correlated to the spent fuel reprocessing as an decreased lifetime increases the number of final (PWR) or start-up (SFR) cores, leading to an increased amount of available spent fuel.

Inventories

Figure 10: Sensitivity of front-end cycle results to each parameter

Total Pu in cycle is increased by having more separated Pu in storage, which in turn is impacted by breeding ratio, SFR introduction, etc. A larger growth rate increases the total fleet size and total material inventory of the fleet, including the Pu in cycle. As the energy growth is exponential, the impact is non-linear and is bigger on the higher side of the sensitivity analysis. The SFR reprocessing parameters also impact the Pu in cycle in a non-linear way as, when a decreased SFR reprocessing results in broken scenarios and in the use of an external source of Pu, it increases the total Pu in the fuel cycle.

The final two tornado diagrams show the impact of input parameters on the total Pu and MA in waste. In both cases, the breeding ratio is the most significant input parameter. As explained previously, the breeding ratio is negatively correlated to the SFR spent fuel flux at reprocessing, which directly impacts the flux of materials being sent to the waste.

Sensitivity table

The sensitivity table provides a way to display the sensitivity analyses results, with all inputs and outputs together on a single The linearity is checked by calculating the coefficient of determination (r 2 ). The sensitivity indicator is retained as a valid parameter when r 2 is greater than 0.9. It must be noticed that each sensitivity indicator is calculated from a small amount of points (three points for most sensitivity indicators), so that the results should be considered as tendencies rather than absolute results.

Compared with the Tornado Diagram that displays two values of sensitivity (qmin and qmax), the sensitivity indicator is more global as it is based on a least squares linear regression upon the three available points (reference, upper sensitive case and lower sensitive case).

On the sensitivity table (Figure 11), a red sensitivity coefficient corresponds to a positive correlation and a blue coefficient to a negative one. A darker colour (red or blue) indicates a higher sensitivity of the output to the parameter. In cases where the relationship is non-linear or otherwise hard to quantify into a single value, the shading is replaced with a "?".

The sensitivity table provides a big picture view of the analysis results. The most important input parameters from a sensitivity perspective can easily be identified by noting the rows with the strongest shading. In accordance with the Tornado diagrams, these are the overall growth rate (energetic production) and the introduction date of the FRs, followed by the rate of FR introduction and the reactor lifetime. Explanations for these impacts have been given in the description of Tornado diagrams.

The change in reactor lifetime impacts material flows without changing the rest of the system.

Additional input parameters with strong impacts on the scenario results are in the areas of breeding gain and reprocessing, where rows have some strongly shaded cells, but also have a number of "?" cells, indicating possible non-linear responses.

Figure 3 :Figure 4 :

 34 Figure 3: Base case scenario -PWR UOX spent fuel flow at reprocessing

Figure 5 :Figure 6 :

 56 Figure 5: Base case scenario -SFR spent fuel flow at reprocessing

Figure 7 :

 7 Figure 7: Base case scenario -Pu inventory in cycle (left) and MA inventory in waste (right)

  1) ref p input parameter's reference value, s p input parameter's value in the sensitivity case ref R / s R output indicator's corresponding to ref p / s p

Table 1 : Identified key parameters and their range of variation

 1

	Parameters	Units	PWR UOX	SFR
	Fissile BU	GWd/tHM	40, 50, 60	100, 136
	Fresh fuel	235 U enrichment	%	4.95 (adjusted with BU)	-
	Equivalent Pu content	%	-	13.8 (adjusted with BU)
	Cycle length	EFPD	410 (adjusted with BU)	340 (adjusted with BU)
	Breeding gain		-	0.75, 0.9, 1, 1.1, 1.25
	Reactor lifetime	y	Infinite, 60, 40	Infinite, 60, 40
	Total nuclear energy demand	TWh/y	430	430 (steady), ,
	Minimum cooling time	y	2, 5, 8	2, 5, 8
	Fabrication time	y	1, 2, 3	1, 2, 3
	Introduction date of FR	y	-	Year 70, 80, 90, 130
	Rate of introduction	y	-	Over 20, 30, 40 years
	First year of reprocessing	y	25, 35, 45	85, 95, 105
	Annual reprocessing capacity	tHM	700, 850, 1000	400, 600, 800
	Losses (U and Pu)	%	0.05, 0.1, 0.2
	Reprocessing priority	-	FIFO to FILO
	Enrichment tail	%	0.15, 0.25, 0.35
	60 -Reference value of the parameter

  table. The correlation of each input (table row) to each output (table column) is shown through colour shading of each cell in the table.A sensitivity indicator (S, see equation (2)) aims at quantifying the variation of an output indicator to the variation of an input parameter, provided the linearity of the output indicator as a function of the input parameter.

	S	=	R p	ref ref	.	p R ∂ ∂	(2)
	ref p input parameter's reference value, p ∂ its variation
	ref R / R ∂ output indicator's corresponding to ref p /	p p ref ∂ +

Figure 11: Representation of the sensitivity study's results with a Sensitivity table

The sensitivity table also reveals more localized high impacts where one or two cells are strongly shaded relative to the rest of a row. For example, the enrichment tails only impacts Unat consumption and SWU requirements, but does not impact the rest of the fuel cycle.

Shifting to the columns, only a few outputs show significant sensitivity to a majority of the inputs, and none are sensitive to all inputs. We can note that outputs that are the most impacted are the storages and inventories.

Concluding remarks

This study aimed at assessing and quantifying the importance of input parameters in scenario analysis. The expert group identified 17 input parameters and evaluated the sensitivity of 22 of the most important scenario outputs to those parameters.

Backend

Backend

Uranium

PWR cycle FR cycle

Reprocessing

Inventory

Front-end Front-end S (r 2 ≥ 0,9) Pu Different scenario codes were used for this study. Where prior benchmarking by this group and others has mainly focused on the common code capabilities, this study used some capabilities only found in some codes. This allowed to include a broader range of cases and to provide fuller information on the sensitivity of the scenario outputs to its inputs.

While the content of this study is most directly applicable to fuel cycle systems analysts, it should also be of value to researchers and developers working on the different parts of the nuclear energy system by providing a systems perspective for their work.


Samuel Bernard
email: bernard@math.univ-lyon1.fr

From Single-Scale to Multiscale Modelling in Biology

Keywords: Multiscale model, master equation, reaction-diffusion equation, hybrid modelling, individual-based modelling, Structured PDE PACS 02.30.Hq, 02.30.Jr, 02.50.Ga Mathematics Subject Classification (2000) 60H10, 35Q84, 60J25, 82C31, 91B69, 35K57

Multiscale modelling in biology is used to take into account the complex interactions between the different organization levels in living systems. We review several models, from the most simple to the complex ones, and discuss their properties from a multiscale point of view.

Introduction

Multiscale modelling first appeared in material engineering problems. Modelers in biology have adapted multiscale techniques to study cancer and other complex biological systems [START_REF] Schnell | Multiscale modeling in biology new insights into cancer illustrate how mathematical tools are enhancing the understanding of life from the smallest scale to the grandest[END_REF][START_REF] Gatenby | A reaction-diffusion model of cancer invasion[END_REF][START_REF] Spencer | Modeling somatic evolution in tumorigenesis[END_REF][START_REF] Zhang | Multiscale agent-based cancer modeling[END_REF][START_REF] Lesart | On the importance of the submicrovascular network in a computational model of tumour growth[END_REF][START_REF] Powathil | Modelling the effects of cell-cycle heterogeneity on the response of a solid tumour to chemotherapy: Biological insights from a hybrid multiscale cellular automaton model[END_REF]. Here we present a few approaches to multiscale modelling in biology. We identify desirable properties such as description at many scales and interactions among and between scales. We proceed stepwise, beginning with classical models in mathematical biology, and introducing progressively richer structures that share more of the features of a multiscale model. We begin with the Lotka-Volterra predator-prey equations [START_REF] Kaplan | Understanding nonlinear dynamics[END_REF]. We then introduce the FKPP reaction-diffusion equation [START_REF] Mendéz | Reaction-transport systems: mesoscopic foundations, fronts, and spatial instabilities[END_REF], followed by the Turing equations for pattern formation [START_REF] Turing | The chemical basis of morphogenesis[END_REF], and the Keller-Segel equations for chemotaxis [START_REF] Keller | Initiation of slime mold aggregation viewed as an instability[END_REF]. We then propose an approach, based on master equations for stochastic processes, to build a multiscale model. We present contemporary approaches for multiscale modelling: individual-based models, hybrid models, and structured PDE models. Examples include multiscale models used to probe cancer growth and treatment strategies [START_REF] Kurbatova | Hybrid model of erythropoiesis and leukemia treatment with cytosine arabinoside[END_REF][START_REF] Ribba | A multiscale mathematical model of cancer, and its use in analyzing irradiation therapies[END_REF][START_REF] Anderson | Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment[END_REF]; models for the cell cycle [START_REF] Bekkal Brikci | An age-and-cyclin-structured cell population model for healthy and tumoral tissues[END_REF][START_REF] Doumic | Analysis of a population model structured by the cells molecular content[END_REF] and cell differentiation [START_REF] Hoffmann | Noise-driven stem cell and progenitor population dynamics[END_REF][START_REF] Bessonov | A multi-agent model describing self-renewal of differentiation effects on the blood cell population[END_REF]; and spatial multi-scale models for tissue and cancer growth [START_REF] Galle | Modeling the effect of deregulated proliferation and apoptosis on the growth dynamics of epithelial cell populations in vitro[END_REF][START_REF] Drasdo | A single-cell-based model of tumor growth in vitro: monolayers and spheroids[END_REF][START_REF] Drasdo | Individual-based models to growth and folding in one-layered tissues: intestinal crypts and early development[END_REF][START_REF] Drasdo | Monte carlo approach to tissue-cell populations[END_REF]. Based on the discussion, we propose two properties any model should possess to be multiscale.

From Single-scale to Multiscale Models

Any working definition of multiscale modelling would classify most mathematical models in biology as non-multiscale, or single-scale. We discuss here the scales of four classical models in biology with increasing complexity.

Predator-Prey Interactions. Lotka-Volterra equations are a system of two nonlinear ordinary differential equations describing the interaction between two species, a predator and a prey [START_REF] Kaplan | Understanding nonlinear dynamics[END_REF]. Their first use in biology date back to 1925 (Lotka) and 1926 (Volterra). The dynamical variables are the populations (numbers) of preys (x) and predators (y). Preys find enough food at all times, and predators feed on prey. The populations evolve according to the equations

x = x(a -by), (1)
y = y(dx -c). (2)
(Here and in the following, time derivatives in ODEs are denoted by the symbol and all coefficient are positive.) At the population scale, equations (1-2) state that the growth rate of the prey population is constant (a) and the death rate of the preys is proportional to the number of predators (by). The death rate of the predators is constant (c) and the growth rate of the predator population is proportional to the prey number (dx). The positive fixed point

(x, y) = ( c d , a b
) is unstable (it is a center). The predator-prey population oscillate around the fixed point along close trajectories (Fig 1). This is an example of a population model, and it involves interactions between objects at a lower scale: individuals. Yet this is clearly not a multiscale model. A multiscale model should be able to provide information at more than one scale, for instance, on population ("is the wolf population declining?") and on the individual themselves ("which wolves find prey?"). The first question can be answered, while the second one cannot. All the wolves are the same, so are all the hares. What is missing to the predator-prey model is the ability to characterize individuals. This example shows that two populations with simple interactions can display oscillations in their number (Fig 1). This is not multiscale modelling because there is only one scale, the population, even though the model rests on assumptions on interactions between individuals. Modeling is about describing interactions between single objects. By considering a large number of these objects, we obtain population equations. It turns out that this link between organization at one scale and the interaction at the lower scale is always present (Fig 2 ).

Reaction-Diffusion Equations. ODEs cannot be multiscale because they do not describe individual objects. Perhaps providing individuals with attributes could help devising a multiscale model. We could endow individual with a spatial position, for instance. Let us have a look at some examples involving a spatial variable x.

The Fisher-Kolmogorov-Petrovskii-Piskunov (FKPP) equation is a reaction-diffusion partial differential equation, introduced by R.A. Fisher in 1937. The equation describes the frequency ρ(t, x) of an advantageous gene in a population located at position x, at time t. Individuals carrying the gene reproduce faster. The growth rate of the frequency ρ is given by a logistic term rρ(1ρ). At the same time, individuals carrying the gene move by diffusion (Dρ xx ). The FKPP reaction-diffusion equation is (Subscripts in PDEs denote partial derivatives.) There are two homogeneous solutions to this equation: an unstable one, ρ(x) ≡ 0, and a stable one, ρ(x) ≡ 1. Equation ( 3) admits traveling wave solutions of the form φ (xct). The wave φ (y) moves at a constant speed c, and connects the steady states 0 and 1, the stable state one "invades" the unstable steady state 0 (Fig 3).

ρ t = Dρ xx + rρ(1 -ρ). (3)
As with the Lotka-Volterra equations, the FKPP equation describes a population. Equivalently, the interaction scale is the individuals carrying or not the advantageous gene. Although we do not see any individual, there are differences in gene frequency depending on the location x. If we look at coarse-grained populations living in different locations given by a partition of the space, like villages dotting the countryside, we can characterize each of them by their gene frequency. Not only can we count the villages, but we can distinguish between them. Therefore, we have satisfied our first multiscale model condition: we have a population model with individual objects (the villages) that can be distinguished. We have two scales of description, the number of villages, and how is each village. Nevertheless, the FKPP as a multiscale model is not satisfying. The number of villages is constant, and the gene frequency does not interact with the number of villages.

Pattern Formation. Alan M. Turing introduced this system in 1952 to show how spatial patterns (spatially heterogeneous solutions) could arise from diffusion of chemical substances, when diffusion was thought to lead to homogeneous solutions. In his paper [START_REF] Turing | The chemical basis of morphogenesis[END_REF], Turing discusses a system of two morphogens regulating each other and diffusing in space. One is an activator and the other is an inhibitor. The activator, u, activates itself and the inhibitor. The inhibitor, v, inhibits itself and the activator. The concentrations of u and v can show spatial instabilities if the diffusion rate of the activator is much smaller than the diffusion rate of the inhibitor. When the activation/inhibition are linear, the equations are

u t = D u u xx + f u u + f v v, (4)
v t = D v v xx + g u u + g v v. (5)
Again, we can coarse-grain the space x to consider, for instance, cells, located at x, characterized by certain amounts of activator and inhibitors. Spatial instabilities, or Turing patterns, offer an example of lasting heterogeneity, by which cells can be distinguished or Chemotaxis. In the multiscale interpretation of the Turing system, cells are static objects. It would be better is they could move and interact, for instance. Chemotaxis is the phenomenon by which cells move according to concentration gradients in their environment. They could be attracted by food, or repelled by poisons. If the cells themselves secrete chemotactic molecules, we can describe the movement of the cells by a model Keller and Segel developed in 1970 [START_REF] Keller | Initiation of slime mold aggregation viewed as an instability[END_REF],

u t = du xx -(uv x ) x , (6)
v t = εv xx + u -av. (7)
The Keller-Segel equations include two scales: cell density (u) and molecule concentration (v). Interaction between cells is based on the chemoattractant concentration they produce (Fig 5). Yet, the two scales are not nested, they live in the same world. Cells and chemoattractant are modelled at the same level, like the predator and the prey in the Lotka-Volterra model. We do not have yet a multiscale model.

Master Equation-Based Models

The Langevin Equation and the Fokker-Planck Equation. To build a multiscale model, we consider y(t) ∈ R the state of a system at time t at a particular scale. If we describe a single cell, y could be mRNA, protein levels, cell mass or gene expression profile. If we describe a population, y could be numbers of quiescent cells, stem cells, differentiated or tumour cells. We assume that y fluctuates around a value given by the nonlinear deterministic system y F(y). We add noise term to get a nonlinear Langevin equation [START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF],

y = f (y) + g(y)ξ (t). (8)
The noise component ξ (t) is usually a white noise, that is, a time-uncorrelated noise with zero mean and finite variance, often chosen Gaussian. The associated Fokker-Planck equation for the probability density of the state y at time t is

P t (y,t) + [ f (y)P(y,t)] y = 1 2 [g(y) 2 P(y,t)] yy . (9)
Fokker-Planck equations and Langevin equations are mathematically equivalent formulations [START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF]. It is also possible to consider Fokker-Planck and Langevin equations in many variables y ∈ R n . As an example, we add a noise term to the Lotka-Volterra system (1-2), which is sensitive to small perturbations. The new Langevin system displays behaviour different from the Lotka-Volterra (Fig 6

The distribution p(y|y ) of jumps in the differentiation space is Gaussian, with a nonlinear, space-dependent variance σ 2 (y ). The rate of jumps is given by R(y), which is assumed to be correlated to the cell proliferation rate r(y). The differentiation domain Γ = R. The associated Fokker-Planck-like equation is found by calculating the first two jump moments [START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF]. The first jump moment, Γ rp(y + r|y)R(y)dr = 0, vanishes, and the second moment, a 2 (y) = Γ r 2 p(y + r|y)R(y)dr = R(y)σ 2 (y). The Fokker-Planck-like equation associated to equation ( 11) is then

n t (y,t) = 1 2 R(y)σ 2 (y)n(y,t) yy + r(y)n(y,t).
Denote by N A (t) the number of cells with y ∈ A, where A ⊆ R is a subset of cells of interest: N A (t) = A n(y,t)dy. We introduce a logistic growth in the following way:

n t (y,t) = 1 2 R(y)σ 2 (y)n(y,t) yy + r(y) 1 - N A (t) K n(y,t). ( 12
)
K is the carrying capacity for population A. Equation ( 12) is not Fokker-Planck equation anymore since it includes a nonlinear population growth term. Cell number is regulated at the cell population scale through the term N A (t), while the single-cell scale defines how cells move in Γ . The associated single-cell Langevin equation is y = R(y)σ (y)ξ (t). The Langevin equation only describes part of equation ( 12):cell movement in Γ . Equation ( 12) has many desirable properties of a multiscale model. It describes a system at two organization scales: the cell (with y ∈ Γ ) and the population (with N A (t)), and the interaction between the two scale through the logistic term. The two scales are nested, N A is formed directly by the density n.

Individual-Based Models

Individual-based modelling (IBM), or agent-based modelling, has been used in computer science, social sciences, ecology and more recently in biology [START_REF] Railsback | Agent-based and individual-based modeling: A practical introduction[END_REF]. Models are composed of many agents who can make decisions, learn and adapt, and interact with other agents and the environment. They are not necessarily multiscale. However, they are flexible enough to allows a multiscale description. Drasdo and colleague have developed a methodology for simulating proliferating cells [START_REF] Drasdo | Monte carlo approach to tissue-cell populations[END_REF]. Cells are modelled as oriented, deformable spheroids that can move, rotate, change shape, grow in volume and divide into two daughter cells. Cells interacts mechanically and biochemically through membrane surface molecules. Internal properties of the cell regulate how fast cells grow and divide, and which type of surface molecules are expressed. Although simulations can be computationally costly, this kind of modeling can reproduce many of the features of a growing tissue. This type of single-cell-based model has had many applications: epithelium [START_REF] Drasdo | Monte carlo approach to tissue-cell populations[END_REF][START_REF] Galle | Modeling the effect of deregulated proliferation and apoptosis on the growth dynamics of epithelial cell populations in vitro[END_REF], tumors [START_REF] Drasdo | A single-cell-based model of tumor growth in vitro: monolayers and spheroids[END_REF][START_REF] Ramis-Conde | Multi-scale modelling of cancer cell intravasation: the role of cadherins in metastasis[END_REF] and more recently liver regeneration [START_REF] Hoehme | Prediction and validation of cell alignment along microvessels as order principle to restore tissue architecture in liver regeneration[END_REF]. A related method has been proposed [START_REF] Newman | Many-body theory of chemotactic cell-cell interactions[END_REF].

Hybrid Continuous-Discrete Models

This is a popular approach to multiscale modelling, in which cells are discrete entities and molecular concentrations are given by continuous equations. The space contains a coarsedgrained lattice. At each node, one or many cells can be present. Each cell is represented individually and is endowed with relevant property: shape, intracellular state, cycling status, mutation, etc. Cells can move, replicate, die and interact with other cells, directly, or via an extracellular medium.

Ribba and colleagues proposed a multiscale model for cancer growth, with the purpose of optimizing therapeutic irradiation protocols [START_REF] Ribba | A multiscale mathematical model of cancer, and its use in analyzing irradiation therapies[END_REF]. The model incorporates gene regulation, cell kinetics, tissue dynamics, macroscopic tumor evolutions and radio-sensitivity dependence on cell cycle phase. Intracellular interactions are modelled by a boolean network. The state of the boolean network determines how cells progress in the cell division cycle and their fate. Cell can either divide, stop in a quiescent phase, or die by apoptosis. Cell cycle progression, arrest or death is monitored by a cell cycle status. Cell are laid on a lattice and are subject to a changing environment, consisting in the local cell density and the oxygen concentration. Cell fate depends on the local environment. On a tissue level, a fluid dynamics model is used to describe cell movement. Radiation therapy affects the molecular network, which in turn affect cell fate, and tumor progression.

Although the model is too complex to reproduce here, we can identify the multiscale features of the model. The model spans three organization scales: genetic/molecular networks, cells and tissues and interactions between scales are modelled explicitly.

Anderson and colleagues developed a multiscale model for tumor morphology and phenotypic evolution, in which phenotypic mutations and selection drive the tumor morphology evolution [START_REF] Anderson | A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion[END_REF][START_REF] Anderson | Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment[END_REF]. Discrete cell are laid on a square lattice. Cells are characterized by a life cycle governing proliferation and death, and by a phenotype with traits describing key properties of cancer cells like propensity to proliferate or cell-cell adhesion. Cell can undergo phenotypic mutations that will affect their ability to proliferate and move. Because resources (space, oxygen) are limited, phenotypic selection operates and defines the dynamics of tumor evolution. Oxygen and extracellular matrix concentrations are modelled with PDEs.

More recently, we have developed a hybrid multiscale model cell proliferation and differentiation [START_REF] Kurbatova | Hybrid model of erythropoiesis and leukemia treatment with cytosine arabinoside[END_REF][START_REF] Bessonov | A multi-agent model describing self-renewal of differentiation effects on the blood cell population[END_REF]. Cells are modelled as discrete objects (IBM). Cell behaviour depends on continuous intracellular and extracellular processes. The model was developed to simulate the evolution of immature blood cells in the bone marrow. In the version proposed by Kurbatova et al. [START_REF] Kurbatova | Hybrid model of erythropoiesis and leukemia treatment with cytosine arabinoside[END_REF], cell differentiation, self-renewal, and apoptosis are determined by an intracellular network, and by extracellular regulation. The resulting multiscale model was used to probe the effect of different chemotherapeutic treatment schedules on leukemic and normal cells.

A feature of hybrid models is their ability to reproduce relevant biological phenomena, including spatial and intercellular heterogeneities. To illustrate the potential of hybrid models, we show how discrete cells with molecular network of circadian clocks can communicate and synchronize their clocks (Fig 8). The circadian clock model oscillate with a period of around 24h only if the cells are properly coupled [START_REF] Bernard | Synchronization-induced rhythmicity of circadian oscillators in the suprachiasmatic nucleus[END_REF].

Structured PDE Models

Population models structured by molecular content are becoming increasingly popular. Structured equations can be derived from Fokker-Planck equations without noise (for instance, equations ( 8) and ( 9) with g = 0). Doumic presented and analysed a model structured by age a and molecular content x [START_REF] Doumic | Analysis of a population model structured by the cells molecular content[END_REF]. The equation describes the evolution of the density n(t, a, x) at time t, of cells aged a with molecular content x. The evolution of the molecular content x, the structure, and the age are given by a system of ODEs, x = F(a, x) and a = 1. Cells are lost with a rate t=0h t=13.25h t=2.75h t=5h t=7.75h t=10.5h t=16h t=18.75h t=21.5h

Fig. 8 Synchronization of circadian clocks in 300 cells induced by cell-cell contact interaction. Cells (small spheres) move freely, but tend to stay in contact with each other, and express a clock marker (blue: low, red: high levels). Details on the intracellular network and coupling can be found in [START_REF] Bernard | Synchronization-induced rhythmicity of circadian oscillators in the suprachiasmatic nucleus[END_REF].

n(t, a = 0, x) = 2 ∞ 0 Γ
b(a, x, y)n(t, a, y)dyda.

This system of equations involves two scales: the molecular scale given by the equation x = F(a, x) and the population scale, with the death term B(a, x) and the boundary condition. This model is a simplified version of a population model structured by cyclins with explicit cell cycle phases [START_REF] Bekkal Brikci | An age-and-cyclin-structured cell population model for healthy and tumoral tissues[END_REF].

Discussion and Conclusion

We have presented here several approaches to multiscale modelling, from master equationbased approaches to individual-based and hybrid models to structured models. Stochastic and multiscale formalisms share many attributes. Depending on the point of view, the same equation can describe a stochastic or a multiscale phenomenon, as with the Fokker-Planck equation. In both multiscale and stochastic models, there are at least two distinct scales.

The choice of the scales is left to the modeler, but many models incorporate scales ranging from molecules to the tissue, with the cell as the fundamental modelling unit (Table 1). One consequence of multiscale modelling is the need to have a detailed but simple description of the intracellular dynamics, with mechanisms to create heterogeneity between cells. Stochastic processes are compatible with the multiscale framework developed here and have been well studied in the context of intracellular networks [START_REF] Kepler | Stochasticity in transcriptional regulation: origins, consequences, and mathematical representations[END_REF][START_REF] Paulsson | Models of stochastic gene expression[END_REF][START_REF] Kaern | Stochasticity in gene expression: from theories to phenotypes[END_REF][START_REF] Perkins | Strategies for cellular decision-making[END_REF][START_REF] Ribeiro | Variability of the distribution of differentiation pathway choices regulated by a multipotent delayed stochastic switch[END_REF]. In a recent review, Byrne and Drasdo discussed the merits of individual-based models and continuum models of cell populations [START_REF] Byrne | Individual-based and continuum models of growing cell populations: a comparison[END_REF]. We could cover only a small part of the multiscale modelling literature. More approaches to multiscale modelling and examples of applications can be found in recent books [START_REF] Chauviere | Cell Mechanics: From Single Scale-Based Models to Multiscale Modeling[END_REF][START_REF] Deisboeck | Multiscale cancer modeling[END_REF][START_REF] Cristini | Multiscale modeling of cancer: an integrated experimental and mathematical modeling approach[END_REF]. Based on the discussions presented here and elsewhere [START_REF]What exactly is multiscale modeling[END_REF], we propose list of properties a computational or mathematical model should possess to be a multiscale model -At least two nested organization scales. We should be able to distinguish attributes of objects at each scale. The lower scale should be imbedded into the higher scale. Coupled equations for cell densities and molecular concentrations like equations ( 6)- [START_REF] Byrne | Individual-based and continuum models of growing cell populations: a comparison[END_REF]. -Interaction between and among scales. Emerging behaviour (higher scale) of interacting particles (lower scale) is not sufficient, the emergent behaviour should interact with the particles themselves.

Fig. 1

 1 Fig. 1 Lotka-Volterra predator-prey system (1), with a = c = 1, b = 0.01 and d = 0.02, for two different initial conditions: x 0 = 20, y 0 = 20 (black) and x 0 = 100, y 0 = 100 (grey). A) Time evolution of preys and predators numbers. B) Solutions shown in panel A lie on distinct, closed trajectories.

Fig. 3

 3 Fig. 3 Solution to the FKPP equation, with r = D = 1, and initial conditions ρ 0 (x) = 1 if x ≤ 0 and 0 otherwise, and Neumann boundary conditions (no flux). A) Raster plot of the solution ρ(t, x), with increasing values from black to white. B) Solution profiles ρ(t, x) in x for successive fixed times t. It can be shown that the solution is "moving" to the right at a constant speed c = 2 √ rD = 2 [25]. The initial condition is emphasized in grey.

Fig. 4

 4 Fig. 4 Solution to the Turing equations, with f u = 0.4, f v = -0.5, g u = 1, g v = -0.5, D u = 0.2, D v = 10, with initial conditions u 0 (x), v 0 (x) = 1 + sin(x), 1 + cos(x) and Neumann boundary conditions (no flux). A) Raster plot of the solution u(t, x), with increasing values from black to white. B) Solution in x at successive, fixed times. The initial condition is emphasized in grey. A B

Fig. 5

 5 Fig. 5 Solution to the Keller-Segel equations (6-7), with d = 0.5, ε = 0.05, a = 0.1, with initial conditions u 0 (x), v 0 (x) = 2.5 1 [-0.9,0.1] + 1 [0.1,0.9] , 0 and Neumann boundary conditions (no flux). A) Raster plot of the solution u(t, x), with increasing values from black to white. White indicates u ≥ 4.0. B) Solution in x at successive, fixed times. The initial condition is emphasized in grey. Solution at t = 2.5 is out of bound.

Fig. 6

 6 Fig. 6 Two realizations of the stochastic process given by the Langevin version of the Lotka-Volterra system: x = x(aby) + eξ (t) and y = -y(cdx) + f ξ (t), for two different initial conditions: x 0 = 20, y 0 = 20 (black) and x 0 = 100, y 0 = 100 (grey). Parameters are are the same as simulations shown in Fig 1, a = c = 1, b = 0.01 and d = 0.02. The noise ξ (t) is a Gaussian white noise with standard deviations e = f = 3. Absorbing boundary conditions were imposed at x = 0 and y = 0. A) Time evolution of preys and predators numbers. In one case (black), the predator goes extinct at around t = 55, leaving the prey population to grow exponentially. B) Solutions shown in panel A wander around. The black preys grow exponentially once the predators are extinct (horizontal trajectory at y = 0).

Fig. 7

 7 Fig. 7 Solution to the multiscale reaction-diffusion equation (12) derived from the Fokker-Planck equation. The domain is Γ = [0, 1], and A= [y N , 1]. Functions are R(y) = r 0 r 1 r(y), r(y) = 4r max y(1y/y N )1 [0,y N ] , σ 2 (y) = s 0 + 4s 1 (y -1/2)2 . Parameter values are r 0 = 0.01, r1 = 0.20, r max = 0.01, s 0 = 0.005, s 1 = 0.02, K = 0.05, y N = 0.8. A) Rate functions. B) Raster plot of the solution n(y,t), with increasing values from black to white. C) N [0,1] (t) (grey) and N A (t) (black). D) Some realizations of Langevin equation (grey) and their mean (black) compared to the mean differentiation value of the solution n(y,t) (dashed line).

  Organization scales in biology. Single-scale model equations are determined by the interactions at the lower scale. For example, to model a molecular network it is necessary to describe the interaction between

	ORGANISATION	INTERACTION	FORMALISM
		WHOLE-BODY MODELLING
			PDE, ODE
	TISSUE/ORGAN	POPULATION DYNAMICS
			PDE, ODE,
			IBM
	CELL	NETWORKS OF NETWORKS
			ODE, BN
	NETWORK	MOLECULAR NETWORKS
			ODE, SDE,
			BN
	MOLECULE

MD

Fig.

2

molecules. From this point of view, a tissue is a collection of cells interacting together; a cell is a collection of networks interacting together; a molecular network is a collection of molecules interacting together. Scales can go lower (atoms/ions interacting in a single molecule) or higher (organisms interaction in an ecosystem). ODE: Ordinary Differential Equations, PDE: Partial Differential Equations, SDE: Stochastic Differential Equations, IBM: Individual-Based Modelling, BN: Boolean Network, MD: Molecular Dynamics.

  ). In particular, population can now go extinct.The Master Equation. Langevin equations can be used to simulate Markov processes. Markov processes are memoryless stochastic processes, and the probability density function of y at time t, P(y,t) can be expressed as a differential equation local in time but nonlocal in space Γ .

	P t (y,t) =	{W (y|y )P(y ,t) -W (y |y)P(y,t)}dy .	(10)
		Γ
	Equation (10) is a master equation. The master equation is a gain-loss equation for the
	probability of each state y. (The Fokker-Planck equation is a special kind of master equation
	used as an approximate description of a Markov process in which jumps are small and the
	nonlinearities are smooth.) As an illustration of the master equation approach, we study a
	model in which cell differentiate by diffusing (i.e. non-directed random movement) into a
	differentiation space y ∈ R [18]. The cell density n(y,t) obeys an equation derived from a
	master equation,
	n t (y,t) =	p(y|y )R(y )n(y ,t)dy -R(y)n(y,t) + r(y)n(y,t).
	Γ

Table 1

 1 Model scales and formalisms.

	Model		Organization Scale
		Mol. Network Cell Kinetics Tissue Dynamics Ext. Space
	Master Eq.	ODE/SDE	Phase	PDE/FP-based	N/A
	Ribba [32]	BN	Phase	PDE	PDE
	Drasdo [15]	N/A	Age	IBM	N/A
	Kurbatova [23]	ODE	Phase	IBM	PDE
	Anderson [3]	N/A	Phen./Phase	CA	PDE
	Bekkal [4], Doumic [11] ODE	Phase/Age	PDE	N/A

ODE: Ordinary Differential Equations, PDE: Partial Differential Equations, Age: Age-structured model, Phase: Phase-structured model, Phen.: Phenotypic space with random mutations, IBM: Individual-Based Modeling, CA: Cellular Automaton, BN: Boolean Network.

B(a, x) and a boundary conditions at a = 0 defines the birth rate of new cells. The structured equation is a birth-death transport equation describing how cells move, die and are born in the age/structure space (a, x). n t + n a + {F(a, x)n} x + B(a, x)n = 0,

Acknowledgements I would like to thank the organizers of the SFBT (Francophone Society for Theoretical Biology) for the opportunity to present an early version of this work. I also thank colleagues from the Inria Dracula Team for insightful discussions, and C. Knibbe for code development used in Fig 8.


Catherine Foley
email: foley@math.mcgill.ca

Samuel Bernard
email: s.bernard@biologie.hu-berlin.de

Michael C Mackey
email: mackey@cnd.mcgill.ca

Cost-effective G-CSF therapy strategies for cyclical neutropenia: mathematical modelling based hypotheses

Keywords: cyclical neutropenia, G-CSF treatment, mathematical modelling, non-standard treatment

Using computer simulations of a previously published mathematical model for the regulation of stem cell and neutrophil production in dogs, we have studied the efficacy of three different treatment protocols for cyclical neutropenia involving granulocyte colony stimulating factor (G-CSF). We used a model based on canine data to infer treatments for humans because of the availability of extensive dog data and the similarity of the dog and human disorders. The first treatment scheme involves alternate day administration of G-CSF. The second triggers G-CSF administration whenever neutrophil levels fall below a predetermined level, and the third one follows a random administration protocol. The computer simulations predict that clinically desirable results can be achieved with all three methods, using far less G-CSF than would be needed with the standard daily treatment. If the results of this modelling are borne out clinically, it would entail a considerable financial savings to patients.

Introduction

Cyclical neutropenia (CN) is a rare haematological disease characterized by oscillations in the circulating neutrophil count. These levels fall from normal to barely detectable levels with a typical period of 19 to 21 days in humans [START_REF] Guerry | Periodic hematopoiesis in human cyclic neutropenia[END_REF][START_REF] Dale | Cyclic neutropenia: A clinical review[END_REF][START_REF] Haurie | Cyclical neutropenia and other periodic hematological diseases: A review of mechanisms and mathematical models[END_REF], even though periods up to 40 days have been observed [START_REF] Haurie | Cyclical neutropenia and other periodic hematological diseases: A review of mechanisms and mathematical models[END_REF]. These oscillations in the neutrophil count are generally accompanied by oscillations around normal levels in other blood cell lineages such as platelets, lymphocytes and reticulocytes [START_REF] Haurie | Cyclical neutropenia and other periodic hematological diseases: A review of mechanisms and mathematical models[END_REF][START_REF] Haurie | Modeling complex neutrophil dynamics in the grey collie[END_REF]. All grey collies [START_REF] Lund | Cyclic neutropenia in grey collie dogs[END_REF] are born with this congenital disease and oscillations with period on the order of 11 to 16 days are observed [START_REF] Haurie | Cyclical neutropenia and other periodic hematological diseases: A review of mechanisms and mathematical models[END_REF](Haurie et al , , 1999b[START_REF] Haurie | Modeling complex neutrophil dynamics in the grey collie[END_REF]. This animal model has provided extensive experimental data to decipher the nature of CN.

Many mathematical models have been proposed to explain the origin of these oscillations as well as to understand the control of neutrophil production in non-pathological cases. For a discussion of previous models that have been developed see [START_REF] Hearn | Cyclical neutropenia and the peripheral control of white blood cell production[END_REF] and [START_REF] Haurie | Cyclical neutropenia and other periodic hematological diseases: A review of mechanisms and mathematical models[END_REF]. In a recent study [START_REF] Bernard | Oscillations in cyclical neutropenia: New evidence based on mathematical modeling[END_REF], it was shown that all of the characteristics of CN in grey collie can be accounted for by an elevated level of cellular death (apoptosis) in the neutrophil precursors. This elevation of apoptosis in CN was first predicted over twenty years ago [START_REF] Mackey | Unified hypothesis of the origin of aplastic anemia and periodic hematopoiesis[END_REF] and has been observed clinically [START_REF] Aprikyan | Impaired survival of bone marrow hematopoietic progenitor cells in cyclic neutropenia[END_REF]. The modelling predicts that decreasing the level of apoptosis in neutrophil precursors can attenuate or eliminate many of the symptoms of CN.

CN is often treated using granulocyte colony stimulating factor (G-CSF) [START_REF] Hammond | Treatment of cyclic neutropenia with granulocyte colony stimulating factor[END_REF] which is known to interfere with apoptosis [START_REF] Koury | Programmed cell death (apoptosis) in hematopoiesis[END_REF][START_REF] Park | Cytokine regulation of apoptosis in hematopoietic precursor cells[END_REF][START_REF] Williams | Haemopoietic colony stimulating factors promote cell survival by suppressing apoptosis[END_REF][START_REF] Williams | Molecular regulation of apoptosis: Genetic controls on cell death[END_REF]. Treatment protocols (Amer. Soc. Clin. Oncol., 1997) typically call for daily subcutaneous injection of G-CSF at 5 µg per kg of body weight, and for a 60 kg adult this currently entails a yearly cost of approximately US$40,000 (CF: Check D.C. Dale reference for other dosage AM.J Hematol. 2003). Clearly it would be of enormous help if the same clinical effects could be achieved with the use of less G-CSF. A few alternative treatment strate-gies have been reported in which various administration schemes have been used [START_REF] Jayabose | Recombinant human granulocyte colony stimulating factor in cyclic neutropenia: Use of a new 3-day-a-week regimen[END_REF][START_REF] Dicato | Reduced G-CSF dosage in cyclic neutropenia gives satisfactory clinical laboratory and economic results[END_REF][START_REF] Danielson | Intermittent g-csf in cyclic neutropenia[END_REF], yet no attempt to use current knowledge of the dynamic behaviour of the haematopoietic system for designing new treatment protocols has apparently been made.

This paper uses the recent model of [START_REF] Bernard | Oscillations in cyclical neutropenia: New evidence based on mathematical modeling[END_REF] to examine, through computer simulations, the predicted effects of different protocols for the administration of G-CSF to grey collies. The goal is to see whether one can find ways of administering G-CSF less frequently than daily and still achieve the same clinical effect currently reached using a daily regimen. In the first section, we describe the model for cyclical neutropenia and its dynamic characteristics. Second, we illustrate the computer predicted results of three different schemes for the administration of G-CSF. The first involves giving G-CSF every other day. The second is a "threshold triggered decision" in which G-CSF is administered whenever the neutrophil count falls below a predetermined value. Finally, the third scheme is a random delivery protocol. The paper concludes with a brief discussion.

The Model

The model used in this paper considers two cell types whose dynamics are regulated by negative feedback loops. These cell types are the haematopoietic stem cells (HSCs) and the mature neutrophils. HSCs are assumed to be self-renewing and multi-potential. That is, HSCs are self-sustained and newly divided cells can either remain in the HSC pool, or be triggered to undergo differentiation into any of the blood cell lineages.

Fig 1 . illustrates the two components of this model: the HSC compartment and the mature neutrophil compartment. The rate at which HSCs differentiate into the neutrophil lineage is assumed to be determined by the blood G-CSF level. As these neutrophil precursors differentiate, their number is amplified through successive divisions.

As shown in Fig 1 ., there are two feedback loops. The first loop, between the mature neutrophil compartment and the HSC compartment, is the rate (F (N )) of HSC differentiation into the neutrophil lineage. F (N ) is mediated by G-CSF levels and operates with a lag of duration τ N that accounts for the time required for neutrophil precursor division and maturation.

The second feedback loop regulates the rate (K(S)) at which HSCs enter the proliferative cycle from the resting (G 0 ) state. It operates with a delay τ S that accounts for the cell cycle duration, that is the length of time required to produce two daughter HSCs from one mother cell. K(S) regulates the level of haematopoietic stem cells, while F (N ) controls the neutrophil count. G-CSF acts on model parameters in various ways. First, it decreases the apoptosis rate of neutrophil precursors and HSCs, leading to an increase of HSC levels and neutrophil amplification [START_REF] Basu | Evaluation of role of g-csf in the production, survival, and release of neutrophils from bone marrow into circulation[END_REF]. Second, it increases the differentiation rate F of HSC into the neutrophil precursor lineage. Third, the HSC cell cycle duration and the neutrophil precursor transit time are decreased. G-CSF is partly cleared through binding to neutrophil receptors. Thus, G-CSF clearance decreases as the number of neutrophils decreases [START_REF] Kearns | Disposition of recombinant human granulocyte colony stimulating factor in children with severe chronic neutropenia[END_REF][START_REF] Terashi | Close association between clearance of recombinant human granulocyte colony-stimulating factor (G-CSF) and G-CSF receptor on neutrophils in cancer patients[END_REF] and the neutrophil count increases when the level of G-CSF is increased [START_REF] Petros | Pharmacokinetics and administration of colony-stimulating factors[END_REF][START_REF] Chatta | Effects of in vivo recombinant methionyl human granulocyte colony stimulating factor on the neutrophil response and peripheral blood colony forming cells in healthy young and elderly adult volunteers[END_REF][START_REF] Price | Effect of recombinant granulocyte colony stimulating factor on neutrophil kinetics in normal young and elderly humans[END_REF]. This type of regulation regulation is a negative feedback mechanism. A similar mechanism is assumed to be in place to control HSC levels, and we take K(S) to be also a negative feedback loop.

In the context of the mathematical modelling, this interaction has been kept as simple as possible to make model-based predictions possible. The equations of the model are presented in Appendix A.

Behaviour of the model

The model predicts that for elevated apoptosis rates, oscillations in both the HSC and neutrophil popula-tions appear [START_REF] Bernard | Oscillations in cyclical neutropenia: New evidence based on mathematical modeling[END_REF]. More interestingly, it was shown that for wide parameter range, the system displays bistability. In this case, depending on the history of neutrophil and HSC levels, the neutrophil count may either oscillate or stabilize to a steady state. This bistability phenomenon suggests a potential neutrophil count or phase-dependent treatment on neutrophil levels. Although there is no experimental evidence that bistability really occurs in CN patients, some experimental records (Haurie et al , 1999b) display different responses to G-CSF treatment applied at different times.

Fig 2 . shows the absolute neutrophil count (ANC) from a grey collie over 712 days. Two phases of continuous G-CSF treatment were applied (day 0 to day 96 and day 296 to day 436). The two phases without treatment show significant differences in the amplitude of oscillation of ANC and periodicities. Spectral analysis using the Lomb periodogram (a generalization of Fourier spectral analysis for unevenly sampled data) has proven to be useful for analysis of periodic haematological diseases [START_REF] Fortin | Periodic chronic myelogenous leukemia: Spectral analysis of blood cell counts and etiological implications[END_REF]Haurie et al , 1999a). We applied the Lomb periodogram technique to the two phases without treatment. Phase 1 (day 96 to 296) has a period of 14.6 days with strong significance p ≤ 10 -5 . Phase 2 also shows some periodicity to the eye (period 15.4 days) but it is not significant. This clear change in the neutrophil dynamics may indicate that the second G-CSF treatment stabilized the neutrophil level, whereas the first phase of treatment did not. Such bistability, if true, would suggest that the neutrophil count can either display periodic fluctuations with large amplitude or stabilize around an average value. In this paper, based on computer simulations we show how the timing of the treatment affects the behaviour of the ANC.

Investigating administration patterning

Daily G-CSF treatment generally results in a marked increase in average neutrophil count and an elevation in the nadir counts, decreasing the duration of severe neutropenia. An absolute neutrophil count (ANC) below 0.5×10 9 cell/L usually defines neutropenia. In investigating different patterning of G-CSF administration with this mathematical model, the number of days of neutropenia in a cycle serve as a benchmark for evaluating efficiency of different treatment strategies.

The mathematical model describes 2 dynamical variables and contains 11 physiological parameters, including 2 time delays. The detailed description of the model is presented in Appendix A. By varying some of the parameters, we can simulate the effects of G-CSF administration in CN. In Appendix B, we explain how we mimicked the administration of G-CSF as well as the numerical methods that have been used.

The effects of non-standard G-CSF delivery protocols

To study the effects of non-standard delivery protocols for G-CSF, we compared a daily administration protocol with alternate day, neutrophil-leveldependent and random administration protocols. For each scheme, computer simulations were performed over 70 days (5 cycles of 14 days). In each case, we computed the solution for every starting day in the cycle (14 runs for each treatment scheme). The average number of neutropenic days as well as the average number of days of treatment are reported in Table I. Daily G-CSF administration is generally associated with an increase of the amplitude and a decrease of the period of the oscillations (Haurie et al (1999a)). However, it has also been reported that in some isolated cases G-CSF therapy abolished significant oscillations (Haurie et al (1999a), [START_REF] Hammond | Treatment of cyclic neutropenia with granulocyte colony stimulating factor[END_REF]).

Computer simulations of the model predict that daily administration of G-CSF could either lead to large amplitude oscillations with reduced period of neutropenia, or stabilization of the ANC above the neutropenic level (average of 1.45×10 9 cell/L, see Fig 3 .). We found that the type of response depends on the neutrophil level on the first day of treatment. This suggests that the starting date of the treatment will affect the behaviour of the neutrophil count. For daily and alternate day therapies, this bistability can be characterized in terms of starting day of the treatment (see Fig 4 .) and both behaviours (stabilization and oscillations) are distinguished in Table I. The treatment would be optimal when the mathematical model predicts a stabilization of the ANC, as shown in Fig 4 . From our numerical simulations, if day 0 represents the day in the cycle at which the ANC nadir occurs, the optimal phase for starting the treatment is between day 3 and day 8 for a cycle length of 14 days. Moreover, the model predicts that starting the treatment during the neutropenic phase (days 0-1 and days 12-13) will lead to large oscillations.

Alternate day delivery

In the every-other-day therapy, the model predicts that the blood neutrophil level when G-CSF administration starts (and hence, the starting day of the treatment) also has an effect on the behaviour of the neutrophils. However, because the treatment is interrupted every two days, the behaviour of the neutrophil count is influenced by the ANC every time G-CSF is injected. Consequently, a more erratic response results, but having the same general behaviour, and we can still distinguish a stabilization and an oscillatory phase (see ANC remains between 0.7 and 2.3 ×10 9 cells/L. The range of ideal starting days is as broad as for daily G-CSF administration, and this scheme also efficiently increases the nadir of the neutrophil level and shortens the period of neutropenia. Therefore, it seems to achieve positive results similar to daily G-CSF. In particular, one can see from Table I that the results for alternate-day treatment during the stabilization phase are better than the results for the oscillating period of the daily scheme. This suggests that if one is able to localize the stabilization phase for a subject, then starting every-other-day treatment in this period would have similarly good results as continuous G-CSF administration in the oscillatory phase. And more importantly, it uses half the amount of G-CSF.

Threshold triggered delivery

To determine the effects of threshold triggering drug delivery, we simulated G-CSF administration whenever the ANC falls below a predetermined level. Results for levels 0.7, 0.9 and 1.2×10 9 cell/L are shown in Table I. The bistability in the parameter space was still present in the results, but there was no obvious distinction between a stabilization and oscillatory phase as a function of starting day of treatment, as for the alternate day and daily schemes (see Fig 6 .). This protocol targets specific days when G-CSF is needed and, as a result, it seems that its efficiency is good with respect to both number of neutropenic days and number of treatment days, especially for high predetermined levels (see for instance level 1.2 ×10 9 cell/L in Table I). Although this scheme presents very in- teresting theoretical results, it is not a practical approach since it would require monitoring of the blood count frequently and hence, involves expensive laboratory work.

Random protocol

The results obtained with the level dependent scheme justify finding a scheme that would give similar results, but that would be more easily implemented. We consider a random protocol in which G-CSF is administered randomly for a predetermined number of days. In particular, we use the number of days of administration obtained in the level dependent method as a first indicator (average of 53.8% of treatment days for level 1.2×10 9 cell/L). However, because the random protocol is obviously not as targeted as the level-dependent scheme, it appears that the number of days of treatment must be higher than for threshold triggered therapy to achieve analogous results. We simulated random administration of G-CSF during 70 days for two cases (50% and 70% of treatment days) (see Fig 7 .). Because the random protocol does not necessarily administrate G-CSF the first day of treatment, we only made simulations for starting day 0 and 5 (in the oscillatory and stabilization phase respectively). In each case, we carried out 12 runs and averaged the results, presented in Table I. The protocol with 50% of treatment days did not give good results. However, the results for the protocol with 70% of treatment days give similar results as the alternate day therapy. Table I gives a summary of all the methods and their results.

Discussion

Cyclical neutropenia has been extensively modelled in dogs due to its interesting dynamics and the availability of extensive data for grey collies. However, there have been no attempt to use its dynamic features to improve existing treatment strategies. Using a mathematical model based on canine data, we have provided a computer simulation study of three treatment protocols that achieve essentially the same results as daily administration of G-CSF but with smaller total G-CSF levels and/or fewer days of administration. The first treatment scheme is the alternate day (or every-other-day) regimen, which increases the nadir of the neutrophil level to clinically acceptable levels. The second proposed treatment protocol, neutrophil level-dependent administration, in its current form would be difficult to implement, but an alternative random treatment scheme could easily be designed to achieve similar results. Mimicking the G-CSF administration pattern that was obtained by computer simulation is a possibility that should be considered. However, because there may be some subtle differences between the dog model and the human disease, clinical trials will be necessary to establish efficacy of the treatments in humans. Moreover, the existence of a region of bistability could be exploited by perturbing the haematopoietic system to stabilize ANC levels. If this were possible, cyclic neutropenia could be triggered into intermittent phases of low level oscillatory, or non-oscillatory, behaviour (c.f. Fig 5 . for an example) and only occasional doses of G-CSF would then be required to maintain this state. This idea of using a "physiological black hole" to trigger or destroy oscillations was originally put forward by [START_REF] Winfree | Time and timelessness in biological clocks[END_REF][START_REF] Winfree | The Geometry of Biological Time[END_REF] in a different context.

Because of the oscillatory nature of cyclical neutropenia and the nonlinear interactions between G-CSF and circulating neutrophils through clearance processes, this disorder is a good candidate for designing dynamic treatment protocols. We have shown, based on a computer simulated model, that phasedependent protocols should be able to achieve equal or better results than standard protocols, while using less G-CSF for fewer days of treatment. Obviously, this result opens up the possibility of decreasing sideeffects as well as increasing the cost-effectiveness of the treatment, thereby improving the quality of life for patients.

The delayed arguments account for the time required for cell division cycles and maturation to occur. The net loss rates are:

L S = S(t)[F (N (t)) + K(S(t))],
(5) L N = αN (t).

(6)

The variable S represents the number of haematopoietic stem cells in the resting (G 0 ) phase. Cells in the resting phase can either enter the proliferative phase at a rate K(S) or differentiate at a rate F (N ) to ultimately give rise to mature neutrophils, represented by the second variable, N . Cells in the HSC proliferative phase undergo apoptosis at a rate γ S and the cell cycle duration is τ S . The fraction of HSC surviving cell division is then exp(-γ S τ S ). Successive divisions amplify cell numbers in the differentiation pathway by a factor A, which is also used to account for cell loss due to apoptosis. After a time τ N , each differentiated HSC gives rise to A mature neutrophils N , which are released into the circulation. It is assumed that mature neutrophils are eliminated at a rate α.

Two feedback loops control the entire process through the neutrophil differentiation rate F (N ) and the HSC proliferation rate K(S). G-CSF mediates neutrophil differentiation through a clearance mechanism. When neutrophil count decreases, C-CSF level goes up due to the lack of clearance and differentiation increases. As a function of the neutrophil count, the differentiation rate is thus monotonically decreasing (here G-CSF is implicitly taken into account), and is described by a decreasing Hill function,

F (N ) = f 0 θ n 1 θ n 1 + N n , (7)
where f 0 is the maximal differentiation rate (i.e. in absence of neutrophil), n is a cooperativity coefficient and θ 1 is the neutrophil concentration at which differentiation is half maximal. Notice that θ 1 depends on endogenous G-CSF production and G-CSF administration. HSC proliferation is assumed to be regulated through a similar mechanism. When the HSC level is depleted, intercellular signals trigger cell replication. These signals become fainter when HSC level increases back to normal levels. The HSC proliferation rate feedback term is thus also described by a decreasing Hill function,

K(S) = k 0 θ s 2 θ s 2 + S s . (8)
The parameter k 0 is the maximal proliferation rate, s is a cooperativity coefficient and θ 2 the HSC concentration at which HSC proliferation is half maximal. See [START_REF] Bernard | Oscillations in cyclical neutropenia: New evidence based on mathematical modeling[END_REF] for a detailed derivation of the feedback functions. For more details on the choice and estimation of those parameters, see [START_REF] Bernard | Oscillations in cyclical neutropenia: New evidence based on mathematical modeling[END_REF]. Numerical simulations were performed using the DDE solver dde23 [START_REF] Shampine | Solving ddes in, applied numerical mathematics[END_REF] for Matlab, which solves delay differential equations with constant delays. Mimicking G-CSF injection in a cyclic neutropenic subject was done by simply changing the corresponding sets of parameters. Five parameters were modified under the influence of G-CSF: the neutrophil amplification A, the HSC apoptosis rate γ S , the half-differentiation level θ 1 , the cell cycle duration τ S and the maturation transit time τ N . We first computed the solution with the CN parameters for a transient period of 200 days. This allowed the system to settle down to the stable periodic solution and gives a good approximation of the behaviour in the CN situation. Then, administration of G-CSF was simulated by solving the system using the set of parameters corresponding to G-CSF. For instance, to simulate the alternate day treatment, we first integrated the DDE system using the CN parameters for one day. Then, we integrated the equations using the G-CSF parameters for the next day, and we kept switching from one parameter set to the other during the desired simulated treatment period. Note that the values of the delays are modified to mimic G-CSF. Because dde23 only allows constant delays, this leads to assumptions on G-CSF kinetics since we could not simulate a gradual decrease in the time delays and the other parameter values when switching from CN to G-CSF. This corresponds to assuming that the effects of G-CSF are immediate when injected, and that those effects last exactly one day before disappearing completely and instantaneously.

  Fig 1: G-CSF mediated neutrophil differentiation.

  Fig 2: (A) Absolute neutrophil count from a grey collie (Oprah, data generously supplied by Prof. David C. Dale). The red bars show periods of continuous G-CSF treatment. Notice the periodicity in region 1 versus region 2. (B) and (C) show the corresponding Lomb periodograms of phases 1 and 2 respectively.

  Fig 3: Stabilization of an oscillatory neutrophil count when starting the daily treatment right after the beginning of the falling phase of oscillation.

Fig 4 :

 4 Fig 4: Effect of starting day of treatment for cycles with a period of 14 days using the daily and alternate day schemes. Ideal starting day of treatment (in the stabilization phase) is between day 3 and day 8.

Fig 5 :Fig 6 :

 56 Fig 5: Simulation of the ANC and HSC count using an alternate day regimen, representing a stabilization (the ANC remains between 0.7 and 2.3 ×10 9 cells/L). Note the low level neutrophil oscillations following the cessation of treatment.

Fig 7 :

 7 Fig 7: Example of the effects of the random protocol where there were 70% of treatment days.

  Summary of the different protocols. Values represent means, and standard deviations are in parentheses. For the level dependent protocol, levels are in ×10 9 cell/L.

  Successive divisions amplify cell number in the differentiation pathway by a factor A, which is also used to account for cell loss due to apoptosis. After a time τN , differentiated cells become mature neutrophils N and are released into the blood. It is assumed that mature neutrophils die at a fixed rate α. Two feedback loops control the entire process through the proliferation rate K(S) and the differentiation rate F (N ).

	used for mimicking the neutrophil production in three				proliferation S K(S)
	situations of interest: in a healthy subject (without CN), in a dog suffering from CN, and under G-CSF				τ S
	treatment.	G	1	S proliferative phase 2 apoptosis γ G HSC	M	2e γ τ	S	G 0 HSC S
							N τ	amplifi-A
		bone marrow	G-CSF mediated differentiation			cation
					S F(N)
		circulation					N
									mature
									neutrophils
					disappearance	α
		Appendix B. Simulating CN and
		G-CSF administration
		Bernard et al (2003) showed that the model can pre-
		dict cyclical neutropenia and simulate G-CSF admin-
		istration by using physiologically relevant values for
		the 11 different parameters. The parameters used for
		numerical simulation are those presented in Bernard
		et al (2003) and they are in agreement with experi-
		mental data. Table II presents the parameter values

S S Fig 8: Model of neutrophil production. The variable S represents the number of haematopoietic stem cells in the resting (G0) phase. Cells in the resting phase can either enter the proliferative phase at a rate K(S) or differentiate at a rate F (N ) to ultimately give rise to mature neutrophils N , the second variable. Cells in the HSC proliferative phase undergo apoptosis at a rate γS and the cell cycle duration is τS.

Table II :

 II Estimated model parameters used for numerical simulations of healthy condition, cyclical neutropenic condition and under G-CSF treatment. Values from[START_REF] Bernard | Oscillations in cyclical neutropenia: New evidence based on mathematical modeling[END_REF].

	Parameter Unit	Health CN	G-CSF
	A	100	380	10	20
	f 0	day -1	0.8	0.8	0.8
	θ 1	10 8 cell/kg 0.36	0.36	0.8
	k 0	day -1	8.0	8.0	8.0
	θ 2	10 6 cell/kg 0.095	0.095 0.095
	n	-	1	1	1
	s	-	2	2	2
	τ N	day	3.5	3.5	3.0
	τ S	day	2.8	2.8	2.6
	γ S	day	0.07	0.07	0.05
	α	day -1	2.4	2.4	2.4

Acknowledgments

This work was supported by MITACS (Canada) and the Natural Sciences and Engineering Research Council (NSERC grant OGP-0036920, Canada). We thank Prof. David Dale and Ms. Caroline Haurie for providing the grey collie data.

Appendix

A.

Mathematical model

The model for neutrophil production consists of a system of two delay differential equations (DDEs). Each equation considers the net production and loss rates (termed P (t) and L(t) respectively) of haematopoietic stem cells and mature neutrophils.

The net production rates are:

P N (t) = AS(t -τ N )F (N (t -τ N )).

(4)