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

Walid Djema, Catherine Bonnet, Hitay Ozbay, Frederic Mazenc

► To cite this version:

Walid Djema, Catherine Bonnet, Hitay Ozbay, Frederic Mazenc. A practical cell density stabilization technique through drug infusions: a simple pathfinding approach. 2023 European Control Conference (ECC 2023), Jun 2023, Bucharest, Romania. 10.23919/ECC57647.2023.10178254. hal-04379210

HAL Id: hal-04379210

https://inria.hal.science/hal-04379210v1

Submitted on 8 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A practical cell density stabilization technique through drug infusions: a simple pathfinding approach

W. Djema

C. Bonnet

H. Özbay

F. Mazenc

**Abstract**— 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.

**Key Words:** Modeling, Switching, Pathfinding.

## 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., [3], [9], [18]). 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 extraconditions- the existence of a unique strictly positive steady state, which is the most important operating point (see, e.g., [4], [22], [19], [2], [3], [1], [20], [13], [9], [8], [11], [17]). 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 these 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], [4], [20]). 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 subsystems<fnref n="1" /> 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

Catherine Bonnet and Frédéric Mazenc are with Université Paris-Saclay, Inria, CNRS, CentraleSupélec, Laboratoire des signaux et systèmes, Gif-sur-Yvette, France. Walid Djema is with Inria, BIOCORE, Centre Inria d'Université Côte d'Azur, France. Hitay Özbay is with Department of Electrical and Electronics Engineering, Bilkent University, Ankara, Turkey.

Emails: walid.djema@inria.fr, catherine.bonnet@inria.fr, hitay@bilkent.edu.tr and frederic.mazenc@inria.fr

<footnote n="1">¹A set of distinct nonlinear systems with distributed delays ([3], [24]), 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 [8], obtained from the design of suitable strict Lyapunov-Krasovskii functionals (LKFs), to characterize the stable steady state of each subsystem.</footnote>
out switching rules that achieve healthy hematopoiesis".
Now we address this issue through simple planning tools
([23], [16]). 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, \dots, 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 growth-
factor dependent. The relationship between the growth-
factor concentrations $e_j, j = 1, \dots, 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 dif-
ferent thresholds on the growth-factor concentrations
$e_j$, the corresponding controlled biological parameters
($K_i, \gamma_i, \beta_i(0), \text{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 parameter-
transitions. This bring us into the fundamental field of
*automated planning and scheduling*, which belongs to
the branch of *artificial intelligence* (AI) [23]. Mainly it
concerns the elaboration and the realization of strategies
and action sequences, in order to set goals and achieve
them (see [16]). This theory is used when the solutions
are complex and must be *discovered* and *optimized* in
multidimensional spaces ([16], [23]). 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 reg-
ulating 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 typi-
cally 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 $\bar{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 neighbor$^2$,
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], [24], [3]). Our approach implicitly
uses classical *planner algorithms*, such as the A-star
algorithm ($A^*$, [15], [21]). Due to space limitation, only
unhealthy hematopoiesis is discussed in this paper.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Example of step-like functions modeling the behavior of the system parameters with respect to their growth-factor concentrations; $\beta(0)$ is a proliferation parameter, $\gamma$ stands for apoptosis, $K$ represents cell differentiation, the self-renewing parameter $L$ is also affected since $L = 1 - K$; see [12].</figcaption></figure>

## 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,

<footnote n="2">²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.</footnote>
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 = \frac{2K\Im}{\mu}\beta(R_e)R_e, \quad \text{where, } \beta(\ell) = \frac{\beta(0)}{1 + b\ell^r},$$

with $\Im > 0$, $\mu > 0$, $\beta(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_{eff}$. 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 $\pm 0.02$ at each time it changes continuously with respect to its growth-factor concentration $\epsilon_4$. We focus on the three main parameters: proliferation $\beta(0)$, apoptosis $\gamma$, 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_e^{(0,0,0)}$, i.e. the activated initial subsystem has the parameters $(\beta(0)^{(0)}, \gamma^{(0)}, K^{(0)})$, then the point denoted$^3 M_e^{(1,0,+1)}$ belongs to the neighborhood of the initial node. However, it is not possible to make a direct transition from $M_e^{(0,0,0)}$ to $M_e^{(0,0,+2)}$, 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 growth-factor dependent parameters has been introduced and

analyzed in [3]. The positive steady state $(R_e, M_e)$ in that case exists if and only if the HSC-compartment satisfies the condition, $(2L \int_0^\tau e^{-\gamma^a} f(a) da - 1) \beta(0) > \delta$, where $f$ is a probability density functions s.t. $\int_0^\tau f(a) da = 1$, (see [3]). 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_{eff}$) 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_{eff}$). 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.

<figure><img src="image_3.png" /><figcaption>Fig. 2. An oscillatory behavior of an unstable node in one subsystem belonging to the overall model. HSC-compartment (blue) parameters: $K = 0.05$, $\gamma = 0.03$, $\tau = 2.81165$, $\delta = 0.85$, $f(a) = \frac{e^a}{e^a - 1}$, and, $\beta(R) = \frac{8}{1+R^3}$. Mature cells (red) parameter: $\mu = 0.025$.</figcaption></figure>

Moreover, because of the high toxicity levels caused by the drugs (e.g. side-effects of chemotherapy), we

<footnote>$^3 M_e^{(+1,0,+1)}$ is a compact form that means that, starting from $M_e^{(0,0,0)}$, the $\beta(0)$ increases by one level, the apoptosis rate remains unchanged, while differentiation increases by one jump.</footnote>

<footnote>$^4$We refer to a subsystem that does not satisfy the exponential stability condition given [8] as "unstable", even if this subsystem can be stable in fact. Indeed, the stability conditions used to determine the candidate nodes are sufficient conditions.</footnote>
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).

<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

## 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 [5], [6]). 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 [15]. 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^*$ [21], for tracking a moving (cell density) reference $M^*(t)$ for all $t \ge 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^* \in \text{Goal}$.

We are focusing on the case where three parameters can be switching: $\gamma$, $\beta(0)$ and $K$ (Fig. 1). Thus, we consider a 3D space denoted xyz, which is in fact the space $\beta(0)-\gamma-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 3D-representation 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 $\beta(0)-\gamma-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 parameter-space 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_{eff}$.

<figure><img src="image_5.png" /><figcaption>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 insignificant (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.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>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_e^{(attribute_1, attribute_2, attribute_3)}$.</figcaption></figure>

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

The concept of heuristics is by definition problem-specific. For instance, in our application, we notice that the nonlinear relationship between the controlled parameters and the steady states ($(\beta(0), \gamma, K) \rightarrow M_e(\beta(0), \gamma, 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_e^{(attribute_\beta, attribute_\gamma, attribute_K)}$ in the effective hematopoietic network, $M_{eff}$, characterized by $G \in Goal$, $G^{(attribute_\beta, attribute_\gamma, attribute_K)}$, is given by,

$$H\left(M_e^{(attribute_\beta, attribute_\gamma, attribute_K)}\right) = \frac{1}{\Delta_\beta} \left| attribute_\beta - attribute_\beta^G \right| + \frac{1}{\Delta_\gamma} \left| attribute_\gamma - attribute_\gamma^G \right| + \frac{1}{\Delta_K} \left| attribute_K - attribute_K^G \right|$$

where $\Delta_\beta$, $\Delta_\gamma$ and $\Delta_K$, are respectively the amplitudes of the one-step (an elementary switching) transition in the characteristic step-like patterns $e_2 \rightarrow \beta(e_2, 0)$, $e_3 \rightarrow \gamma(e_3)$, and, $e_4 \rightarrow K(e_4)$. For instance in Fig. 1 we note that $\Delta_\beta = 0.3$, $\Delta_\gamma = 0.05$ and $\Delta_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).

<figure><img src="image_7.png" /></figure>

<table><caption>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.</caption><tr><td>The inadmissible node ● can be:<ul><li>an unstable steady state</li><li>non-existing positive steady state</li><li>a point outside the Limit_area</li></ul></td><td>The heuristic has to be:<ul><li>Admissible: It must always be lower or equal to the actual path.</li><li>Consistent: H(N<sub>1</sub>) ≤ cost(N<sub>1</sub>, N<sub>2</sub>) + H(N<sub>2</sub>)</li></ul></td></tr><tr><td colspan="2">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 fully-determined through their distance from the prescribed reference M* and satisfy the Confidence_area requirement.</td></tr></table>

## III. NUMERICAL ILLUSTRATIONS

We give a numerical example that illustrates the overall stabilization technique based on Algorithm 1. The fixed parameters of the system are as follow. For HSCs: $\tau = 1.25$ (cell-cycle duration), $\delta = 0.14$ (death rate), $f(a) = \frac{me^{ma}}{e^{m\tau}-1}$ (mitosis function) where $m = 10$, $\beta(e_2, R) = \frac{\beta(e_2, 0)}{1+bR^n}$ (re-introduction function from resting to proliferating phases [18]), with $b = 1$, $n = 4$. For mature cells the degradation rate is $\mu = 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 ([3]) is not verified. Furthermore, only 59 subsystems have a positive steady state that satisfies the exponential stability condition ([8]). We set $\text{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. $\text{length}(M_{eff}) = 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 [10]), the drug has only a limited effect on $\beta(0)$. Thus, we consider without loss of generality that $\beta$ can vary only between 0.95 and 2.15. It follows that the number of nodes in the network is reduced to 22 ($\text{length}(M_{eff}) = 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 $\beta(0)$-values to those given in Table I, then a third candidate $M_e$ exists in Goal, for the triplet: $\beta(0) = 2.45$, $\gamma = 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.

<figure><img src="image_8.png" /><figcaption>Fig. 7. The hematopoietic network corresponding to Table I.</figcaption></figure>

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:

<figure><img src="image_9.png" /></figure>

The corresponding therapeutic strategy, i.e. the model-parameters 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], [10], [13], [2]. 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 [5], [6], [7]), to better characterize the optimal times at which the drug should be administrated.

## REFERENCES

[1] Adimy, M., Bourfia, Y., Hbid, M.L. and Marquet, C., 2016. Age-structured model of hematopoiesis dynamics with growth factor-dependent coefficients. Elect. J. Diff. Eq., p.140.

[2] Adimy, M. and Crauste, F., 2007. Modelling and asymptotic stability of a growth factor-dependent stem cells dynamics model with distributed delay. Discrete and Continuous Dynamical Systems-Series B, 8(1), pp.19-38.

[3] Adimy, M., Crauste, F., ElAbdllaoui, A., 2008. Discrete-maturity structured model of cell differentiation with applications to acute myelogenous leukemia. Journal of Biological Systems, 16(03), pp.395-424.

[4] Bélair, J., Mackey, M.C. and Mahaffy, J.M., 1995. Age-structured and two-delay models for erythropoiesis. Mathematical biosciences, 128(1-2), pp.317-346.

[5] Djema, W., Bayen, T. and Bernard, O., 2022. Optimal Darwinian Selection of Microorganisms with Internal Storage. Processes, 10(3), p.461.

[6] Djema, W., Giraldi, L., Maslovskaya, S. and Bernard, O., 2021. Turnpike features in optimal selection of species represented by quota models. Automatica, 132, p.109804.

[7] Djema, W., Bernard, O. and Giraldi, L., 2020. Separating two species of microalgae in photobioreactors in minimal time. Journal of Process Control, 87, pp.120-129.

[8] Djema, W., Mazenc, F. and Bonnet, C., 2017. Stability analysis and robustness results for a nonlinear system with distributed delays describing hematopoiesis. Systems & Control Letters, 102, pp.93-101.

[9] Djema, W., Mazenc, F., Bonnet, C., Clairambault, J., Hirsch, P. and Delhommeau, F., 2016. Stability of a delay system coupled to a differential-difference system describing the coexistence of ordinary and mutated hematopoietic stem cells. Conference on Decision and Control (CDC), pp. 561-566.

[10] Djema, W., Bonnet, C., Mazenc, F., Clairambault, J., Fridman, E., Hirsch, P., and Delhommeau, F. (2018). Control in dormancy or eradication of cancer stem cells: Mathematical modeling and stability issues. Journal of Theoretical Biology, 449, 103-123.
<table><caption>TABLE I</caption><thead><tr><th>Nodes</th><th>β(0)</th><th>γ</th><th>K</th><th>M<sub>e</sub></th><th>M<sub>e</sub><sup>(0)</sup></th><th>M* ± 0.5</th></tr></thead><tbody><tr><td>1</td><td>0.95</td><td>0.40</td><td>0.06</td><td>0.81808465</td><td>✓</td><td></td></tr><tr><td>2</td><td>0.95</td><td>0.40</td><td>0.08</td><td>0.98614342</td><td></td><td></td></tr><tr><td>3</td><td>0.95</td><td>0.35</td><td>0.12</td><td>1.70355101</td><td></td><td></td></tr><tr><td>4</td><td>0.95</td><td>0.35</td><td>0.14</td><td>1.35798982</td><td></td><td></td></tr><tr><td>5</td><td>0.95</td><td>0.30</td><td>0.16</td><td>2.45099792</td><td></td><td></td></tr><tr><td>6</td><td>0.95</td><td>0.30</td><td>0.18</td><td>2.47511923</td><td></td><td></td></tr><tr><td>7</td><td>0.95</td><td>0.25</td><td>0.20</td><td>3.24792752</td><td></td><td></td></tr><tr><td>8</td><td>0.95</td><td>0.25</td><td>0.22</td><td>3.41423221</td><td></td><td></td></tr><tr><td>9</td><td>1.25</td><td>0.45</td><td>0.06</td><td>0.89302665</td><td></td><td></td></tr><tr><td>10</td><td>1.25</td><td>0.40</td><td>0.10</td><td>1.77626758</td><td></td><td></td></tr><tr><td>11</td><td>1.25</td><td>0.35</td><td>0.14</td><td>2.66862097</td><td></td><td></td></tr><tr><td>12</td><td>1.25</td><td>0.35</td><td>0.16</td><td>2.76695795</td><td></td><td></td></tr><tr><td>13</td><td>1.25</td><td>0.30</td><td>0.20</td><td>3.93743744</td><td></td><td></td></tr><tr><td>14</td><td>1.55</td><td>0.45</td><td>0.06</td><td>1.26462747</td><td></td><td></td></tr><tr><td>15</td><td>1.55</td><td>0.45</td><td>0.08</td><td>1.46594557</td><td></td><td></td></tr><tr><td>16</td><td>1.55</td><td>0.40</td><td>0.12</td><td>2.65308368</td><td></td><td></td></tr><tr><td>17</td><td>1.55</td><td>0.35</td><td>0.18</td><td>3.65098593</td><td></td><td></td></tr><tr><td>18</td><td>1.85</td><td>0.45</td><td>0.08</td><td>2.00965592</td><td></td><td></td></tr><tr><td>19</td><td>1.85</td><td>0.40</td><td>0.14</td><td>3.52953970</td><td></td><td></td></tr><tr><td>20</td><td>1.85</td><td>0.30</td><td>0.24</td><td>4.48839554</td><td>✓</td><td></td></tr><tr><td>21</td><td>2.15</td><td>0.45</td><td>0.10</td><td>2.74132806</td><td></td><td></td></tr><tr><td>22</td><td>2.15</td><td>0.40</td><td>0.14</td><td>4.33497547</td><td>✓</td><td></td></tr></tbody></table>

The 22 nodes constituting the vector $M_{eff}$, with their $β(0)-γ-K$ coordinates, i.e. attributes.

<table><caption>TABLE II</caption><thead><tr><th>Nodes</th><th>β(0)</th><th>γ</th><th>K</th><th>M<sub>e</sub></th><th>M<sub>e</sub><sup>(0)</sup></th><th>Goal</th><th>Heuristic</th></tr></thead><tbody><tr><td>1</td><td>0.95</td><td>0.40</td><td>0.06</td><td>0.81808465</td><td>✓</td><td></td><td>8</td></tr><tr><td>2</td><td>0.95</td><td>0.40</td><td>0.08</td><td>0.98614342</td><td></td><td></td><td>7</td></tr><tr><td>3</td><td>0.95</td><td>0.35</td><td>0.12</td><td>1.70355101</td><td></td><td></td><td>6</td></tr><tr><td>4</td><td>0.95</td><td>0.35</td><td>0.14</td><td>1.35798982</td><td></td><td></td><td>5</td></tr><tr><td>5</td><td>0.95</td><td>0.30</td><td>0.16</td><td>2.45099792</td><td></td><td></td><td>7</td></tr><tr><td>6</td><td>0.95</td><td>0.30</td><td>0.18</td><td>2.47511923</td><td></td><td></td><td>8</td></tr><tr><td>7</td><td>0.95</td><td>0.25</td><td>0.20</td><td>3.24792752</td><td></td><td></td><td>10</td></tr><tr><td>8</td><td>0.95</td><td>0.25</td><td>0.22</td><td>3.41423221</td><td></td><td></td><td>11</td></tr><tr><td>9</td><td>1.25</td><td>0.45</td><td>0.06</td><td>0.89302665</td><td></td><td></td><td>8</td></tr><tr><td>10</td><td>1.25</td><td>0.40</td><td>0.10</td><td>1.77626758</td><td></td><td></td><td>5</td></tr><tr><td>11</td><td>1.25</td><td>0.35</td><td>0.14</td><td>2.66862097</td><td></td><td></td><td>4</td></tr><tr><td>12</td><td>1.25</td><td>0.35</td><td>0.16</td><td>2.76695795</td><td></td><td></td><td>5</td></tr><tr><td>13</td><td>1.25</td><td>0.30</td><td>0.20</td><td>3.93743744</td><td></td><td></td><td>8</td></tr><tr><td>14</td><td>1.55</td><td>0.45</td><td>0.06</td><td>1.26462747</td><td></td><td></td><td>7</td></tr><tr><td>15</td><td>1.55</td><td>0.45</td><td>0.08</td><td>1.46594557</td><td></td><td></td><td>6</td></tr><tr><td>16</td><td>1.55</td><td>0.40</td><td>0.12</td><td>2.65308368</td><td></td><td></td><td>3</td></tr><tr><td>17</td><td>1.55</td><td>0.35</td><td>0.18</td><td>3.65098593</td><td></td><td></td><td>5</td></tr><tr><td>18</td><td>1.85</td><td>0.45</td><td>0.08</td><td>2.00965592</td><td></td><td></td><td>5</td></tr><tr><td>19</td><td>1.85</td><td>0.40</td><td>0.14</td><td>3.52953970</td><td></td><td></td><td>1</td></tr><tr><td>20</td><td>1.85</td><td>0.30</td><td>0.24</td><td>4.48839554</td><td></td><td></td><td>8</td></tr><tr><td>21</td><td>2.15</td><td>0.45</td><td>0.10</td><td>2.74132806</td><td></td><td></td><td>3</td></tr><tr><td>22</td><td>2.15</td><td>0.40</td><td>0.14</td><td>4.33497547</td><td>✓</td><td></td><td>0</td></tr></tbody></table>

Nodes, switching parameters, and heuristic measures. The heuristic is computed by considering that $M^*$ is the 22nd node.

<figure><img src="image_10.png" /><figcaption>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).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>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.</figcaption></figure>

[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 Conference, 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.
**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 $\epsilon_2 \to \beta(\epsilon_2, 0)$, $\epsilon_3 \to \beta(\epsilon_3)$, $\epsilon_4 \to \beta(\epsilon_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;

10. if Goal = $\emptyset$;

11. return Failure;

12. else

13. if $M_e^{(0)}$ is not a stable node then

14. Find in the neighborhood of $M_e^{(0)}$ the cheapest path leading to a stable node $M_s$;

15. Save the therapeutic action leading from $M_e^{(0)}$ to $M_s$ then consider that $M_e^{(0)} \leftarrow M_s$;

16. else

17. Continue

**Real-time tracking and processing:**

18. for $i \leftarrow 1$ to length(Goal) do

19. Initialize the Open_list with $M_e^{(0)}$;

20. Initialize the Closed_list with unstable points, stable points $\notin$ Limit_area, and points in Goal that are different from Goal($i$); The cost of switching from $M_e^{(0)}$ to itself is set to $G(M_e^{(0)}, M_e^{(0)}) = 0$;

21. Compute the (Manhattan-like distance) heuristic between Goal($i$) and each stable steady state which is inside the limit_area and outside the confidence_area;

22. while Open_list $\neq \emptyset$ do

23. Current_target takes the point that has the lowest $F$ value ($F = H + G$);

24. Discover the Successors which are the neighbors of Current_target;

25. Ignore the neighbors in Successors that already belong to Closed_list;

26. Add the first-time discovered nodes in Successors to the Open_list;

27. Pop off Current_target from Open_list and add it to Closed_list;

28. for $j \leftarrow 1$ to length(Successors) do

29. if $G(\text{Current\_target}) + G(\text{Current\_target}, \text{Successors}(j)) \geq G(\text{Successors}(j))$ then;

30. Continue

31. else

32. (Current_target $\leftarrow$ Successors($j$)). Cost(node)=$Cost$;

33. Therapeutic_strategy$_i$ = [$M_e^{(0,0,0)}$, Best_neighbors$_1^i$, Best_neighbors$_2^i$, ..., Goal($i$)];

34. Therapeutic_strategy$_i$ $\leftarrow$ the best strategy between the *i*-strategies (Therapeutic_strategy$_i$);

35. return Therapeutic_strategy;

36. Time out during the treatment ($T_{\text{treatment}}$ days);

37. Monitoring the patient's conditions after treatment was completed;

38. if Post-therapeutic analysis reveals that $M_e^{\text{post\_therapy}}$ is stable and respects the Confidence_area then

39. return 'The therapy was successful';

40. else

41. Consider the previous failure and return to Initialization

[13] Fridman, E., Bonnet, C., Mazenc, F. and Djema, W., 2016. *Stability of the cell dynamics in Acute Myeloid Leukemia*. Systems & Control Letters, 88, pp.91-100.

[14] Foley, C. and Mackey, M.C., 2009. *Dynamic hematological disease: a review*. J. Mathematical Biology, 58, pp.285-322.

[15] Hart, P.E., Nilsson, N.J. and Raphael, B., 1968. *A formal basis for the heuristic determination of minimum cost paths*. IEEE trans. on Systems Sc. & Cybernetics, 4(2), pp.100-107.

[16] Kanal, L. and Kumar, V. eds., 2012. *Search in artificial intelligence*. Springer Science & Business Media.

[17] Kareva, I., 2016. *Primary and metastatic tumor dormancy as a result of population heterogeneity*. Biology direct, 11(1), p.37.

[18] Mackey, M.C., 1978. *Unified hypothesis for the origin of aplastic anemia and periodic hematopoiesis*. Blood, 51(5), p.941.

[19] Mahaffy, J.M., Bélair, J. and Mackey, M.C., 1998. *Hematopoietic model with moving boundary condition and state dependent delay: applications in erythropoiesis*. Journal of theoretical biology, 190(2), pp.135-146.

[20] Marciniak-Czochra, A., Stiehl, T., Ho, A.D., Jäger, W., Wagner, W., 2009. *Modeling of asymmetric cell division in hematopoietic stem cells—regulation of self-renewal is essential for efficient repopulation*. Stem cells and development, 18(3), pp.377-386.

[21] Stentz, A., 1994. *Optimal and efficient path planning for partially-known environments*. IEEE International Conference on Robotics and Automation (pp. 3310-3317).

[22] Stiehl, T. and Marciniak-Czochra, A., 2011. *Characterization of stem cells using mathematical models of multistage cell lineages*. Math. & Computer Mod., 53(7-8), pp.1505-1517.

[23] Russell, S., Norvig, P. and Intelligence, A., 1995. *A modern approach. Artificial Intelligence*. Prentice-Hall, Egnlewood Cliffs, 25(27), pp.79-80.

[24] Özbay, H., Bonnet, C., Benjelloun, H. and Clairambault, J., 2012. *Stability analysis of cell dynamics in leukemia. Mathematical Modelling of Natural Phenomena*, 7(1), pp.203-234.