W Djema

Catherine Bonnet
email: catherine.bonnet@inria.fr

H Özbay

Frédéric Mazenc
email: frederic.mazenc@inria.fr

A practical cell density stabilization technique through drug infusions: a simple pathfinding approach

Keywords: Modeling, Switching, Pathfinding

We consider a nonlinear system with distributed delays modeling cell population dynamics, where the parameters depend on growth-factor concentrations. A change in one of the growth factor concentrations may lead to a switch in the corresponding model parameter. Our first objective is to achieve a network representation of the switching system involving nodes and edges. Each node stands for a full-fledged nonlinear system with distributed delays where the parameters are constant. For each node, a stable positive steady state may exist. In the network framework, a change in the growth-factor concentration is interpreted as a transition from one node to another. The objective is then to determine the best switching signal steering the biological parameters over time, making the overall dynamic system moving from one operating mode to another, until reaching a desired stable state. Our method provides a (sub)optimal therapeutic strategy, guiding the density of cells from an abnormal state towards a healthy one, through multiple drug infusions. The drug sequence is deduced from the optimal switching signal provided by a classical pathfinding algorithm, associated with the network representation.

I. Problem Statement

There are still many complex open issues related to the stabilization of cell population densities in healthy and unhealthy situations, as discussed throughout this tutorial session. In this talk, we suggest a new approach to exploit some theoretical results on dynamic systems describing cell population dynamics in order to achieve cell density regulation in hematopoiesis. More precisely, we use known stability conditions of meaningful equilibria in a classical model to design a new framework in which we address the issue of cell density regulation through multiple discrete drug infusions. In the literature, a large part of the models describing cell population dynamics involve fixed parameters representing cell differentiation, self-renewing, proliferation, apoptosis rates, etc. (see, e.g., [START_REF] Adimy | Discretematurity structured model of cell differentiation with applications to acute myelogenous leukemia[END_REF], [START_REF] Djema | Stability of a delay system coupled to a differential-difference system describing the coexistence of ordinary and mutated hematopoietic stem cells[END_REF], [START_REF] Mackey | Unified hypothesis for the origin of aplastic anemia and periodic hematopoiesis[END_REF]). All these parameters depend in fact on some growth factors that regulate the overall biological processes. However, taking into account their dynamics will considerably complicate the resulting models, as well as their mathematical analysis. Earlier models usually admit -under extraconditionsthe existence of a unique strictly positive steady state, which is the most important operating point (see, e.g., [START_REF] Bélair | Agestructured and two-delay models for erythropoiesis[END_REF], [START_REF] Stiehl | Characterization of stem cells using mathematical models of multistage cell lineages[END_REF], [START_REF] Mahaffy | Hematopoietic model with moving boundary condition and state dependent delay: applications in erythropoiesis[END_REF], [START_REF] Adimy | Modelling and asymptotic stability of a growth factor-dependent stem cells dynamics model with distributed delay[END_REF], [START_REF] Adimy | Discretematurity structured model of cell differentiation with applications to acute myelogenous leukemia[END_REF], [START_REF] Adimy | Age-structured model of hematopoiesis dynamics with growth factor-dependent coefficients[END_REF], [START_REF] Marciniak-Czochra | Modeling of asymmetric cell division in hematopoietic stem cells-regulation of self-renewal is essential for efficient repopulation[END_REF], [START_REF] Fridman | Stability of the cell dynamics in Acute Myeloid Leukemia[END_REF], [START_REF] Djema | Stability of a delay system coupled to a differential-difference system describing the coexistence of ordinary and mutated hematopoietic stem cells[END_REF], [START_REF] Djema | Stability analysis and robustness results for a nonlinear system with distributed delays describing hematopoiesis[END_REF], [11], [START_REF] Kareva | Primary and metastatic tumor dormancy as a result of population heterogeneity[END_REF]). Many evidences support the hypothesis under which the hematopoietic system admits multiple operating modes. For instance, the body adapts its blood count (through growth-factor secretion) to face some seasonal allergies, or when dealing with asthma, eczema, and infections. In theses cases, the eosinophil normal count moves from one value to another, depending on body requirements. Moreover, unhealthy hematopoiesis which can be caused by many blood disorders (anemia, cyclic neutropenia, overproliferating malignancies, etc.), may return to normality (i.e., reaching a healthy reference for mature cells denoted M * in the sequel) through drug infusions.

In [12], we discussed a model where almost all the biological parameters are growth-factor dependent and which allows the existence of multiple operating modes within an unified framework. This model involves characteristic patterns (step-like functions) that describe the effect of growth-factor concentrations on the biological functionalities (differentiation, apoptosis, etc.). The resulting system is mainly based on the time-scaling heterogeneity that exists between small hormone-like molecules and cell proliferation ([12], [START_REF] Bélair | Agestructured and two-delay models for erythropoiesis[END_REF], [START_REF] Marciniak-Czochra | Modeling of asymmetric cell division in hematopoietic stem cells-regulation of self-renewal is essential for efficient repopulation[END_REF]). This system is used to describe the hematopoietic system, from stem cells to mature cells [12]. We suggested that the overall system is a possibly switching one, where growth-factor concentrations use an event-triggered control on the involved biological parameters. This approach appears to be efficient in describing both the unhealthy (we assume that drugs act as growth factors) and the healthy biological processes. In particular, the model gives the possibility to existence of several strictly positive steady states. It is a complex system composed by several subsystems1 and switching parameters. We have left in [12] an open problem that has been formulated as follows: "After laying the foundation of this switching model, future work will focus on the role of mature cells in the control exerted on growth factors, in order to find out switching rules that achieve healthy hematopoiesis". Now we address this issue through simple planning tools ( [START_REF] Russell | A modern approach[END_REF], [START_REF]Search in artificial intelligence[END_REF]). For that, we firstly introduce a suitable formulation of the stabilization problem of mature blood cell density in a network framework. The objective is to determine how drugs may manage to steer the overall hematopoietic system, through admissible intermediate states/subsystems, to a healthy desirable operating state. Due to page limitations the reader is referred to [12] for a description of the mathematical model and notation used in the present paper. We focus on the model given in Fig. 1 of [12] that is based on earlier works using coupled PDE-systems of McKendrick-type. It governs the dynamics of gradually immature cell subpopulations residing inside the bone marrow, R i , i = 1, . . . , n, n > 1, together with one type of mature cells M that are active in the bloodstream. Almost all the biological parameters and functions, involved in the age-structured PDEs, as well as their associated boundary conditions, are growthfactor dependent. The relationship between the growthfactor concentrations e j , j = 1, . . . , 5, and the values of their corresponding parameters are approximated by some step-like functions, as described in [12]. This leads to a triggered-event operating mode, such that at different thresholds on the growth-factor concentrations e j , the corresponding controlled biological parameters (K i , γ i , β i (0), etc.) jump from one level to another. Thus, by changing the value of a specific growth-factor concentration (e.g. through drug infusion) the value of its corresponding functional parameter switches (see Fig. 1). The model of interest is then described by the switching system in Table 3 of [12]. In the numerical example given in [12], the switching signal is assumed to be known. The crucial issue that we address in this paper is about the determination of how growth-factor concentrations evolve with respect to time in order to drive the events leading to a suitable parameter-switching sequence. The objective is to steer the model trajectories and reach a target, which is formulated in terms of the total density of mature cells (M * ). We seek a systematic strategy to pursue, in order to provide the suitable switching signal managing the optimal succession of parametertransitions. This bring us into the fundamental field of automated planning and scheduling, which belongs to the branch of artificial intelligence (AI) [START_REF] Russell | A modern approach[END_REF]. Mainly it concerns the elaboration and the realization of strategies and action sequences, in order to set goals and achieve them (see [START_REF]Search in artificial intelligence[END_REF]). This theory is used when the solutions are complex and must be discovered and optimized in multidimensional spaces ( [START_REF]Search in artificial intelligence[END_REF], [START_REF] Russell | A modern approach[END_REF]). Using this approach, we want to specify in our context how the switching occurs between appropriate subsystems candidates in order to achieve a final goal (close to M * ). This task depends on whether subsystems defined for given set of parameters represent unhealthy/healthy hematopoitic systems that have stable positive steady states. The regulating strategy operates under strict conditions, some of them stem from the model (mathematical constraints), while others are introduced to take into account realistic specifications (biological constraints). We consider typically the case where we need to move the total density of a type of mature cells from an initial value M (0) to a new value close to the nominal required density M * in a normal body. Moving from M (0) to M * may require several drug infusions, i.e. the model parameters may take several discrete values until reaching desired combination that produces a stable steady state M close to M * . Note that, among a large number of possible subsystems candidates in the path from M (0) to M * , only a few of them are admissible. This is because many times, when switching from a subsystem to a neighbor2 , the strictly positive steady state of the latter subsystem may be mathematically nonexistent, or it may exist but it turns out later that it is unstable (e.g. Fig. 2), or it can be stable but biologically irrelevant. Thus the question here is: how can we choose among all possibilities the best path (optimal succession of admissible states) to move from M (0) to M * , under practical constraints and criteria. This issue is formulated and solved in an algorithmic and practical way, that uses the stability properties of the overall system ([12], [START_REF] Özbay | Stability analysis of cell dynamics in leukemia[END_REF], [START_REF] Adimy | Discretematurity structured model of cell differentiation with applications to acute myelogenous leukemia[END_REF]). Our approach implicitly uses classical planner algorithms, such as the A-star algorithm (A * , [START_REF] Hart | A formal basis for the heuristic determination of minimum cost paths[END_REF], [START_REF] Stentz | Optimal and efficient path planning for partially-known environments[END_REF]). Due to space limitation, only unhealthy hematopoiesis is discussed in this paper.

In the scale-unit of growth-factor concentration (e 2 ) β(0)

II. Stabilization through drug infusions

Here we set some basic rules (non-unique) that are taken into account in order to establish the optimal discrete drug infusion strategy.

A. Set of rules and practical considerations

We focus on the case where the model involves one immature compartment (HSCs, in the bone marrow) and one mature-cell compartment (e.g. white blood cells, in the bloodstream). We denote by "stable node" or "stable subsystem" a subsystem belonging to the set of all subsystems having a stable positive equilibrium point. This equilibrium point is (R e , M e ), where R e > 0 stands for HSCs and M e > 0 for mature cells. We can refer to the positive steady point using only M e , since M e and R e have similar existence and stability properties (see Remark 1). In fact, we note from Table 3

in [12] that, M e = 2Kℑ µ β(R e )R e , where, β(ℓ) = β(0) 1 + bℓ r ,
with ℑ > 0, µ > 0, β(0) > 0, b > 0 and r > 0. This is convenient in practice because the target is usually formulated with respect to the active mature cells in the bloodstream M * . Thus, "stable node" in the network is fully-defined via its associated M e . The domain of all the possible existing nodes is defined as "Limit_area" that contains only the biologically significant values of M e (e.g. the biologically irrelevant values are excluded). Then, the set of all effective stable M e in a given network is denoted by M ef f . Now, we state the first fundamental rule in our network.

Rule 1: The neighbors of a given node are those which can be directly reached from it through one transition, i.e. which requires no more than one switching in each possibly switching model parameter. Fig. 1 shows the form of the characteristic patterns as step-like functions. For instance, K is switching by levels of ±0.02 at each time it changes continuously with respect to its growth-factor concentration e 4 . We focus on the three main parameters: proliferation β(0), apoptosis γ, and differentiation K (it follows that the self-renewal parameter L = 1-K is also switching). Rule 1 allows one elementary move in each possibly switching parameter. For instance, if the initial node is M (0,0,0) e , i.e. the activated initial subsystem has the parameters β(0) (0) , γ (0) , K (0) , then the point denoted 3 M (+1,0,+1) e belongs to the neighborhood of the initial node. However, it is not possible to make a direct transition from M (0,0,0) e to M (0,0,+2) e , since this action requires a double jump in the K-value. In practice, Rule 1 avoids giving more than one dose of the same drug simultaneously. Indeed, the elementary switching-step is defined as the maximum tolerable dose corresponding to one molecule per infusion (giving more can be toxic). However, Rule 1 allows the use of distinct drugs simultaneously if they are targeting different parameters (a combined-targeted therapy).

Remark 1: The fundamental characteristics of each node (i.e. each subsystem) are: i) Existence of M e : determining if the positive steady state M e of mature cells exists or not for each possible subsystem. For instance, the time-delay version of the model without growthfactor dependent parameters has been introduced and 3 M (+1,0,+1) e is a compact form that means that, starting from M (0,0,0) e , the β(0) increases by one level, the apoptosis rate remains unchanged, while differentiation increases by one jump.

analyzed in [START_REF] Adimy | Discretematurity structured model of cell differentiation with applications to acute myelogenous leukemia[END_REF]. The positive steady state (R e , M e ) in that case exists if and only if the HSC-compartment satisfies the condition, 2L τ 0 e -γa f (a)da -1 β(0) > δ, where f is a probability density functions s.t. τ 0 f (a)da = 1, (see [START_REF] Adimy | Discretematurity structured model of cell differentiation with applications to acute myelogenous leukemia[END_REF]). ii) Stability of M e : we check if the positive steady state (R e , M e ) is stable 4 . We mention that, even if the model and its switching parameters are slightly modified, we only require the existence and the stability properties ((i)-(ii)) of the new model in order to apply the planning algorithms that we will use for stabilization (the technique remains valid).

Next, in therapeutic actions we distinguish between two situations: the initial state of the system can be stable, but evolves around an abnormal blood cell count (e.g. a severe low blood count in a type of mature cells); and the initial system can be unstable and shows an oscillatory behavior as in cyclic blood disorders (e.g. chronic cancer or cyclical neutropenia). Thus, when constructing the effective network regrouping only the allowed nodes (M ef f ) we use the following rule.

Rule 2: When the initial system starts from an unstable node, the first therapeutic step is to drive it to the most quickly reachable stable node. For instance, we note in Fig. 3 that in a general construction the optimal pathway (or even the unique one) towards M * could be the one that goes through some unstable nodes. In order to avoid drowning in some details that lengthen the network description, we limit ourselves to the application of Rule 2. Thus, we focus on networks formed only by stable and biologically meaningful nodes (M ef f ). This restriction is also motivated by the fact that often oscillations in blood counts (Fig. 2) may emerge from unstable nodes, exceeding the tolerable limit of blood cells (denoted Limit_area), which is prohibited. Moreover, because of the high toxicity levels caused by the drugs (e.g. side-effects of chemotherapy), we give priority to any therapeutic strategy that requires minimum quantity of drugs at each time-infusion. Thus, we state the following practical rule.

Rule 3: i) The cost of a one-step transition that requires targeting the three biological parameters in one dose, at the same time, is 6u. ii) The cost of a one-step transition that requires targeting two biological parameters in one dose, at the same time, is 3u. iii) The cost of a one-step transition that requires targeting only one biological parameter in one dose is 1u. In Rule 3, the u indicates a normalized unit associated to the cost in our context. The lowest cost (1u) is associated to the simplest transition from one node to another, that requires exactly one switching in only one of the parameters. Clearly, the choice of the tuning weighting constants {1, 3, 6} is not unique. Rules 1-3 defines the outlines of a reasonable framework, by taking into account some practical/biological constraints. These rules can be easily extended (e.g. to include different costs for different drugs, to include different toxicity levels, price or a shortage of some molecules, etc). Fig. 3. Stabilization through drug infusions: an example involving nodes and transitions with their costs (Rule 3). In this illustration, the dark nodes are unstable, while blue ones are stable. Due to space limitation, we consider networks that involve only stable nodes.

M (0)

B. Network representation, planning tools

Generally speaking, there are two main ways to solve the optimal planning problems: dynamic programming and greedy approaches. The former one is based on the idea that the whole problem should be divided in several subproblems, then by combining the individual solutions we can determine the optimal solution of the overall problem (see, e.g., the direct optimization methods used in [START_REF] Djema | Optimal Darwinian Selection of Microorganisms with Internal Storage[END_REF], [START_REF] Djema | Turnpike features in optimal selection of species represented by quota models[END_REF]). On the other hand, the greedy paradigm is an algorithmic and heuristic approach making the locally optimal choice at each stage, by taking one greedy choice after another until finding an optimal (global) solution. Our approach will be based on a variant of the A * search algorithm, which belongs to the family of greedy approaches. The conditions of consistency, admissibility and optimality of A * , that allow it to find the optimal path, are discussed for instance in [START_REF] Hart | A formal basis for the heuristic determination of minimum cost paths[END_REF]. Here we briefly comment the choice of the heuristic function related to planning/scheduling and regulating strategies in our biological context. Other extensions are beyond the scope of the present paper; e.g. application of more recent algorithms, like D * [START_REF] Stentz | Optimal and efficient path planning for partially-known environments[END_REF], for tracking a moving (cell density) reference M * (t) for all t ≥ 0 in healthy hematopoiesis.

Definition 1: At a node M e , H(M e ) is a heuristic that estimates the cost of the cheapest path from the node M e to the considered objective node M * ∈ Goal.

We are focusing on the case where three parameters can be switching: γ, β(0) and K (Fig. 1). Thus, we consider a 3D space denoted xyz, which is in fact the space β(0)-γ-K where each axis admits a set of finite discrete values defined according to the evolution patterns of the biological parameters. The resulting set of nodes in the 3D space is illustrated in Fig. 4, where the M (0) node (abbreviation of M (0,0,0) ) is represented in purple color. Here, each node is placed near its possible neighbors (26 neighbors for a center point, if all of them satisfy the fundamental requirements, Remark 1). However, this 3Drepresentation is not intuitive and does not bring the complete information about the features of the network. It may even be misleading, because some nodes may for instance be in the center of the 3D-construction but their corresponding M e -value is outside the Limit_area. This results from the dependency between the value of M e and the biological parameters of the model, including the switching ones. We prefer a more convenient representation as a conceptual 2D-construction, based directly on the resulting values of M e , as sketched in Fig. 5. More precisely, we set a more intuitive representation where the flat space incorporates all the admissible nodes and excludes the points that are outside the Limit_area. The Confidence_area contains the so-called Goal-nodes, which are the possible objectives (very close to the prescribed reference M * ). Therefore, we can note that the 3D-hematopoietic network represented according to its switching parameters, in the β(0)-γ-K-space, can be transformed into a representation where the objectives are more apparent, related to the total density of mature cells. The nodes in the 2D representation are the stable M e points, the Limit_area defined for the mature cell density, the therapy objectives defined by the prescribed M * density and its Confidence_area (which allow to determine the Goal set). All these elements are defined according to the distance separating the current nodes from the therapy-objective M * . To illustrate a typical situation that can be encountered in the hematopoietic network system, we can consider the case where two nodes belong to the Goal set, i.e. both of them are close to M * (in terms of distance, that quantifies if the total cell-density M is sufficiently close to the target value M * ) and satisfy the Confidence_area requirement (e.g. G 1 and G 2 ), and we highlight that actually G 1 and G 2 can be distant from one another in the parameterspace representation. Such a situation is specific to our application, since in many other applications the points belonging to the Confidence_area are generally also connected -in the vicinity-or neighborhood of each other.

In the sequel, we use the representation in 2D-projection as in the Fig. 5 to illustrate the objectives, however for the implementation of the heuristic function (using A *like algorithms) the 3D-representation as in Fig. 4 is fundamental. Indeed, when it comes to programming the algorithm providing the optimal path, we notice that on one hand the control inputs are the discrete drug infusions (their costs and transitions are defined in the 3D-space of parameters), while on the other hand the Goal set, its Limit_area and effectives nodes are defined with respect to the density of mature cells (as in Fig. 5). Consequently, the pseudocode describing the search techniques exploits both frameworks. Firstly, it defines the parameter-coordinates of all the nodes that can be in Goal, related to the total density of mature cells. Then, these coordinates are used to compute costs and heuristic functions of each node in the network M ef f . Fig. 4. A selected part of the hematopoietic network in the 3D space. The connections between the nodes are not represented. Notice that, in fact, some of the nodes may not exist if the corresponding triplet does not satisfy one of the fundamental requirements. Other nodes can be unstable or stable but insignifiant (outside the Limit_area). However, all this information cannot be deduced directly from the 3D representation, so we are looking for a better representation (a suitable 2D projection), as in Fig. 5.

M(0) M * Fig. 5.
The 2D representation of the network, with some illustrative nodes and connections between them. The costs of the transitions are not mentioned. The nodes outside the Limit_area are not taking into account when searching the optimal path to the objective. The optimal path is the cheapest possible path linking the purple node (initial point) to one of the nodes satisfying the Confidence_area requirements (Goal). Blue nodes are defined by the values of the positive steady states of mature cells and their corresponding triplet of parameters, M (attribute 1 ,attribute 2 ,attribute 3 ) e .

C. What can be the heuristic function in our network?

The concept of heuristics is by definition problemspecific. For instance, in our application, we notice that the nonlinear relationship between the controlled parameters and the steady states ((β(0), γ, K) → M e (β(0), γ, K)), along with the consideration of Rules 1-3, imply that even if the progress of M e points towards the objective M * is measured by the difference between the cell total densities at M * and every node M e in the network, it is not guaranteed that the closest points to M * (in terms of blood cell count) has necessarily the cheapest cost to reach the Goal set (Fig 6). In the unhealthy case, the heuristic function associated with an arbitrary node M (attribute β ,attributeγ ,attribute K ) e in the effective hematopoietic network, M ef f , characterized by

G ∈ Goal, G (attribute G β ,attribute G γ ,attribute G K )
, is given by,

H M attribute β ,attributeγ ,attribute K e = 1 ∆ β attribute β -attribute G β + 1 ∆γ attributeγ -attribute G γ + 1 ∆ K attribute K -attribute G K
where ∆ β , ∆ γ and ∆ K , are respectively the amplitudes of the one-step (an elementary switching) transition in the characteristic step-like patterns e 2 → β (e 2 , 0), e 3 → γ (e 3 ), and, e 4 → K (e 4 ). For instance in Fig. 1 we note that ∆ β = 0.3, ∆ γ = 0.05 and ∆ K = 0.02. The pseudocode that describes our pathfinding technique is given in Algorithm 1 (see, Real-time tracking and processing, this is an adaptation of classical pathfinding A * -type algorithms to our specific model).

G

The heuristic does not depend on the distance between a node and the final destination G. However, the objective points G (belonging to the set Goal) are fullydetermined through their distance from the prescribed reference M* and satisfy the Confidence_area requirement.

S

The inadmissible node can be: • an unstable steady state • non-existing positive steady state • a point outside the Limit_area

The heuristic = cheapest switching through hematopoietic subsystems

The actual switching through stable mature cell nodes of the hematopoetic subsystems.

The heuristic has to be: 1+bR n (re-introduction function from resting to proliferating phases [START_REF] Mackey | Unified hypothesis for the origin of aplastic anemia and periodic hematopoiesis[END_REF]), with b = 1, n = 4. For mature cells the degradation rate is µ = 0.05. Next, using all the possible combinations of the switching parameters given by Figure 1, we define 1000 possible hematopoietic subsystem constituting the overall system. However, for half of them the condition of existence of the positive steady state ( [START_REF] Adimy | Discretematurity structured model of cell differentiation with applications to acute myelogenous leukemia[END_REF]) is not verified. Furthermore, only 59 subsystems have a positive steady state that satisfies the exponential stability condition ( [START_REF] Djema | Stability analysis and robustness results for a nonlinear system with distributed delays describing hematopoiesis[END_REF]). We set Limit_area = [0, 5] (normalized value for the total density of mature cells), the points M e > 5 are excluded since they are not relevant biologically (too high). It follows that the resulting network contains 51 hematopoietic subsystems, i.e. length (M ef f ) = 51. The number of possibilities can be reduced by adding extra-constraints limiting drug toxicity. Moreover, we suppose that due to the presence of epigenetic mutations promoting overproliferation (see Fig. 1 in [START_REF] Djema | Control in dormancy or eradication of cancer stem cells: Mathematical modeling and stability issues[END_REF]), the drug has only a limited effect on β(0). Thus, we consider without loss of generality that β can vary only between 0.95 and 2.15. It follows that the number of nodes in the network is reduced to 22 (length (M ef f ) = 22). It is worth mentioning that the full-system involving 51 nodes can be addressed exactly in the same way.

Remark 2: The nodes 20 and 22 satisfy the Confidence_area requirement and thereby are potential therapy objectives. If we do not limit β(0)-values to those given in Table I, then a third candidate M e exists in Goal, for the triplet: β(0) = 2.45, γ = 0.4, K = 0.16.

According to Table (I), we set the network system representation (Fig. 7) corresponding to our overall hematopoietic system. In the Confidence_area, only the 22nd node is reachable. According to its attributes (i.e. coordinates) the heuristic is computed in Table II. I.

Using Table II and the costs of the edges in Figure 7, we easily determine the optimal path from the initial node to the objective taking into account practical consideration as detailed in the pseudocode (Algorithm 1). The strategy giving the best nodes sequence in this case is given by: The corresponding therapeutic strategy, i.e. the modelparameters evolution based on the switching signal, is the one presented in Fig. 8. Next, we apply the optimal therapeutic strategy to the time-delay system in [12]. The drug infusions are administered starting from t s1 = 50 (a.u), then t s2 = 57, t s3 = 64, t s4 = 71, t s5 = 78, t s6 = 95, and the therapy ends at t s7 = 102. The resulting model trajectories (total densities of HSCs and mature cells) are illustrated in Fig. 9.

IV. Conclusion

Control theory provided strong tools to determine sufficient conditions for the stability of strictly positive steady state of the nonlinear system with distributed delays modelling healthy and unhealthy hematopoiesis [12], [START_REF] Djema | Control in dormancy or eradication of cancer stem cells: Mathematical modeling and stability issues[END_REF], [START_REF] Fridman | Stability of the cell dynamics in Acute Myeloid Leukemia[END_REF], [START_REF] Adimy | Modelling and asymptotic stability of a growth factor-dependent stem cells dynamics model with distributed delay[END_REF]. In particular, techniques based on the design of suitable Lyapunov-Krasovskii functionals allowed us to estimate the size of the basin of attraction of the positive equilibria. The resulting stability conditions are used in this work for stabilization purposes of a more complex dynamic model where the biological functions are growth-factor or drug dependent [12]. Our objective is to determine an admissible optimal control strategy for drug infusions based on classical pathfinding algorithms, that steer the model trajectories of the full-switching system towards a desired state-target (Goal set). Future work will focus on the use of more formal methods, e.g. based on the Pontryagin's Maximum Principle (as in [START_REF] Djema | Optimal Darwinian Selection of Microorganisms with Internal Storage[END_REF], [START_REF] Djema | Turnpike features in optimal selection of species represented by quota models[END_REF], [START_REF] Djema | Separating two species of microalgae in photobioreactors in minimal time[END_REF]), to better characterize the optimal times at which the drug should be administrated. Fig. 8. The desired therapeutic strategy achieving optimal stabilization of the studied hematopoietic system, towards the favorable target M * . Drugs must be infused at each switching-time exactly as it is given in this figure, so that the mature-cell converge towards the desired M * (as in Fig. 9). Fig. 9. Trajectories of the nonlinear system with distributed delays, when applying the therapeutic strategy in Fig. 8, where drug infusions are spread over intervals of one week.

[11] Djema, W., Bonnet, C., Clairambault, J., Mazenc, F., Hirsch, P., Delhommeau, F., 2017. Analysis of a model of dormancy in cancer as a state of coexistence between tumor and healthy stem cells. American Control Boston, USA, p.5135.

[12] Djema, W., Özbay, H., Bonnet, C., Fridman, E., Mazenc, F. and Clairambault, J., 2017. Analysis of Blood Cell Production under Growth Factors Switching. IFAC World Congress, Toulouse, IFAC-PapersOnLine, 50(1), pp.13312-13317.

3Fig. 1 .

 1 Fig.1. Example of step-like functions modeling the behavior of the system parameters with respect to their growth-factor concentrations; β(0) is a proliferation parameter, γ stands for apoptosis, K represents cell differentiation, the self-renewing parameter L is also affected since L = 1 -K; see[12].

Fig. 2 .

 2 Fig. 2. An oscillatory behavior of an unstable node in one subsystem belonging to the overall model. HSC-compartment (blue) parameters: K = 0.05, γ = 0.03, τ = 2.81165, δ = 0.85, f (a) = e a e τ -1 , and, β(R) = 8 1+R 3 . Mature cells (red) parameter: µ = 0.025.

Fig. 6 .

 6 Fig. 6. In this representation, the visually closest points represent the cheapest costs, without having to consider the value of the mature cell steady states.

Fig. 7 .

 7 Fig. 7. The hematopoetic network corresponding to TableI.

  TABLE IThe 22 nodes constituing the vector M ef f , with their β(0)-γ-K coordinates, i.e. attributes. TABLE II Nodes, switching parameters, and heuristic measures. The heuristic is computed by considering that M * is the 22nd node. Drug infusions are spreaded over time (days)

	Nodes	β(0)	γ	K	Me	M	(0) e	M * ± 0.5	Nodes	β(0)	γ	K	Me	M	(0) e	Goal	Heuristic
	1	0.95	0.40	0.06	0.81808465	✓		1	0.95	0.40	0.06	0.81808465	✓		8
	2	0.95	0.40	0.08	0.98614342				2	0.95	0.40	0.08	0.98614342				7
	3	0.95	0.35	0.12	1.70355101				3	0.95	0.35	0.12	1.70355101				6
	4	0.95	0.35	0.14	1.35798982				4	0.95	0.35	0.14	1.35798982				5
	5	0.95	0.30	0.16	2.45099792				5	0.95	0.30	0.16	2.45099792				7
	6	0.95	0.30	0.18	2.47511923				6	0.95	0.30	0.18	2.47511923				8
	7	0.95	0.25	0.20	3.24792752				7	0.95	0.25	0.20	3.24792752				10
	8	0.95	0.25	0.22	3.41423221				8	0.95	0.25	0.22	3.41423221				11
	9	1.25	0.45	0.06	0.89302665				9	1.25	0.45	0.06	0.89302665				8
	10	1.25	0.40	0.10	1.77626758				10	1.25	0.40	0.10	1.77626758				5
	11	1.25	0.35	0.14	2.66862097				11	1.25	0.35	0.14	2.66862097				4
	12	1.25	0.35	0.16	2.76695795				12	1.25	0.35	0.16	2.76695795				5
	13	1.25	0.30	0.20	3.93743744				13	1.25	0.30	0.20	3.93743744				8
	14	1.55	0.45	0.06	1.26462747				14	1.55	0.45	0.06	1.26462747				7
	15	1.55	0.45	0.08	1.46594557				15	1.55	0.45	0.08	1.46594557				6
	16	1.55	0.40	0.12	2.65308368				16	1.55	0.40	0.12	2.65308368				3
	17	1.55	0.35	0.18	3.65098593				17	1.55	0.35	0.18	3.65098593				5
	18	1.85	0.45	0.08	2.00965592				18	1.85	0.45	0.08	2.00965592				5
	19	1.85	0.40	0.14	3.52953970				19	1.85	0.40	0.14	3.52953970				1
	20	1.85	0.30	0.24	4.48839554			✓	20	1.85	0.30	0.24	4.48839554				8
	21	2.15	0.45	0.10	2.74132806				21	2.15	0.45	0.10	2.74132806				3
	22	2.15	0.40	0.14	4.33497547			✓	22	2.15	0.40	0.14	4.33497547			✓	0
									Therapeutic strategy
				2.15
				1.85
			β(0)	1.55
				1.25
				0.95
					Infusion 1	Infusion 2	Infusion 3	Infusion 4	Infusion 5	Infusion 6		Infusion 7
				0.5
				0.45
			γ
				0.4
				0.35
					Infusion 1	Infusion 2	Infusion 3	Infusion 4	Infusion 5	Infusion 6		Infusion 7
				0.16
				0.14
				0.12
			K	0.1
				0.08
				0.06
				0.04
					Infusion 1	Infusion 2	Infusion 3	Infusion 4	Infusion 5	Infusion 6		Infusion 7

A set of distinct nonlinear systems with distributed delays ([START_REF] Adimy | Discretematurity structured model of cell differentiation with applications to acute myelogenous leukemia[END_REF],[START_REF] Özbay | Stability analysis of cell dynamics in leukemia[END_REF]), each one may have an exponentially stable positive steady state under extra conditions on its associated biological parameters. We use the exponential stability results of[START_REF] Djema | Stability analysis and robustness results for a nonlinear system with distributed delays describing hematopoiesis[END_REF], obtained from the design of suitable strict Lyapunov-Krasovskii functionals (LKFs), to characterize the stable steady state of each subsystem.

two subsystems are in the vicinity of each other, if there exists a switching action to pass directly from one subsystem to the other one, and vice versa.

We refer to a subsystem that does not satisfy the exponential stability condition given[START_REF] Djema | Stability analysis and robustness results for a nonlinear system with distributed delays describing hematopoiesis[END_REF] as "unstable", even if this subsystem can be stable in fact. Indeed, the stability conditions used to determine the candidate nodes are sufficient conditions.

Algorithm 1

The general outline of the procedure for the search of an optimal stabilization strategy through multiple drug infusions.

Require: Initial model parameters and explicit dependence patterns between the controlled parameters and their growth-factors concentrations. The reference M * and the time-span between blood tests. Ensure: The optimized therapeutic strategy to stabilize the HSCs and mature cells total densities at a desired final steady state.

Initialization: system parameters, characteristic patterns; 1. The hematopoietic system structure. The set of fixed parameters and functions of the model; 2. The prescribed reference M * , the prescribed Confidence_area, and, Limit_area; 3. The evolution profiles of the controlled functions e 2 → β(e 2 , 0), e 3 → β(e 3 ), e 4 → β(e 4 ); Network construction: 4. From the first medical testings, identify the parameters of the initially activated subsystem; 5. Identify all the hematopoietic subsystems for all the possible parameter combinations; 6. Identify the subsystems that have strictly positive steady states inside the Limit_area; 7. Identify the subclass of systems that have stable strictly positive steady states; 8. Set up the network connections between all the subsystems; 9. Determine the Goal set using the prescribed reference M * and Confidence_area;