
Citation: Nikolai D. Botkin, Varvara L. Turova, Andrey E. Kovtanyuk, Irina N. Sidorenko, Renée Lampe. Extended model of impaired cerebral autoregulation in preterm infants: Heuristic feedback control[J]. Mathematical Biosciences and Engineering, 2019, 16(4): 2334-2352. doi: 10.3934/mbe.2019117
[1] | Peter Rashkov . Modeling repellent-based interventions for control of vector-borne diseases with constraints on extent and duration. Mathematical Biosciences and Engineering, 2022, 19(4): 4038-4061. doi: 10.3934/mbe.2022185 |
[2] | Felicia Maria G. Magpantay, Xingfu Zou . Wave fronts in neuronal fields with nonlocal post-synaptic axonal connections and delayed nonlocal feedback connections. Mathematical Biosciences and Engineering, 2010, 7(2): 421-442. doi: 10.3934/mbe.2010.7.421 |
[3] | Wing-Cheong Lo, Ching-Shan Chou, Kimberly K. Gokoffski, Frederic Y.-M. Wan, Arthur D. Lander, Anne L. Calof, Qing Nie . Feedback regulation in multistage cell lineages. Mathematical Biosciences and Engineering, 2009, 6(1): 59-82. doi: 10.3934/mbe.2009.6.59 |
[4] | Yuqian Mei, Debao Guan, Xinyu Tong, Qian Liu, Mingcheng Hu, Guangxin Chen, Caijuan Li . Association of cerebral infarction with vertebral arterial fenestration using non-Newtonian hemodynamic evaluation. Mathematical Biosciences and Engineering, 2022, 19(7): 7076-7090. doi: 10.3934/mbe.2022334 |
[5] | Yoon-gu Hwang, Hee-Dae Kwon, Jeehyun Lee . Feedback control problem of an SIR epidemic model based on the Hamilton-Jacobi-Bellman equation. Mathematical Biosciences and Engineering, 2020, 17(3): 2284-2301. doi: 10.3934/mbe.2020121 |
[6] | Yannick Lutz, Rosa Daschner, Lorena Krames, Axel Loewe, Giorgio Cattaneo, Stephan Meckel, Olaf Dössel . Modeling selective therapeutic hypothermia in case of acute ischemic stroke using a 1D hemodynamics model and a simplified brain geometry. Mathematical Biosciences and Engineering, 2020, 17(2): 1147-1167. doi: 10.3934/mbe.2020060 |
[7] | Xingjia Li, Jinan Gu, Zedong Huang, Chen Ji, Shixi Tang . Hierarchical multiloop MPC scheme for robot manipulators with nonlinear disturbance observer. Mathematical Biosciences and Engineering, 2022, 19(12): 12601-12616. doi: 10.3934/mbe.2022588 |
[8] | Sue Wang, Jing Li, Saleem Riaz, Haider Zaman, Pengfei Hao, Yiwen Luo, Alsharef Mohammad, Ahmad Aziz Al-Ahmadi, NasimUllah . Duplex PD inertial damping control paradigm for active power decoupling of grid-tied virtual synchronous generator. Mathematical Biosciences and Engineering, 2022, 19(12): 12031-12057. doi: 10.3934/mbe.2022560 |
[9] | Scott R. Pope, Laura M. Ellwein, Cheryl L. Zapata, Vera Novak, C. T. Kelley, Mette S. Olufsen . Estimation and identification of parameters in a lumped cerebrovascular model. Mathematical Biosciences and Engineering, 2009, 6(1): 93-115. doi: 10.3934/mbe.2009.6.93 |
[10] | Kala Agbo Bidi . Feedback stabilization and observer design for sterile insect technique models. Mathematical Biosciences and Engineering, 2024, 21(6): 6263-6288. doi: 10.3934/mbe.2024274 |
Cerebral autoregulation is an important option of the brain vascular system to provide a stable CBF under variations of MAP [1]. This ensures continuous oxygen supply to the brain tissue. Cerebral autoregulation is impaired in preterm infants, i.e., the dependence of CBF on MAP is almost linear, and the autoregulatory plateau usually exists for a very narrow range of MAP [2]. Additionally, the slope of the autoregulatory curve strongly depends on the arterial partial CO2 level (pCO2) [3]. Increased pCO2 values (e.g. permissive hypercapnia), being often observed in ventilated infants, lead to progressively impaired autoregulation. Moreover, variations in intracranial/cerebral venous pressure (ICP/CVP) may cause additional fluctuations of CBF. All this may be damaging for unstable cerebral blood vessels, especially for those of the germinal matrix (GM) [4,5]. The germinal matrix, a highly vascularized site of origin for neuronal and glial cells, disappears by the 32nd week of gestation. Thus, the risk of rupture of unstable GM-vessels is especially high in very preterm infants. Another risk concerns a raised CVP, which can reduce the perfusion pressure defined as the difference between MAP and CVP. This may result in a decrease of CBF with the threat of cerebral hypoperfusion [6,7].
Mathematical modeling of cerebral autoregulation is the topic of many publications, see e.g. [8,9,10,11,12,13,14]. A comprehensive survey on models of CBF autoregulation can be found in the monograph [15]. The purpose of the current paper is to model impaired cerebral autoregulation in premature newborns and to develop a feedback control that prevents large fluctuations of CBF caused by variations of MAP, pCO2, and ICP/CVP. Based on the polynomial autoregulation models proposed by the authors in [14,16], an extended model of impaired cerebral autoregulation, coupled to a model of cerebral blood vessel network from [17], is suggested. The new feature of this model, compared with [14] and [16], consists in accounting for the effect of pCO2 and CVP on CBF. According to the literature (cf. [3,18]), such an effect is very important in ventilated preterm infants. Moreover, a new heuristic feedback control maintaining CBF within a safety range in the presence of unpredictable variations in MAP, pCO2, and ICP/CVP is constructed. Viability theory [19] and differential game approach [20] are used to prove the reliability of this control. It should be noted that the controller designed gives much better results compared with that proposed in [16].
The paper is organized as follows.
Section 2 describes a model of cerebral autoregulation assuming a polynomial dependence of vessel radii on the mean arterial blood pressure. The model is coupled with a hierarchical model of blood vessel network from [17] and adjusted to experimental data for preterm infants using a polynomial fitting.
In Section 3, the autoregulation model from Sect. 2 is modified by adding a mechanism of impaired autoregulation. Variations in MAP, pCO2, and CVP are considered as main factors influencing CBF.
In Section 4, a heuristic feedback controller based on the discrepancy between ideal and impaired autoregulation is proposed. The control action can be interpreted as injection of medicaments causing, with some delay, the dilation or shrinkage of vessels. In Subsection 4.1, a conflict control problem with appropriate state constraints is introduced, and a state-based feedback controller is designed in Subsection 4.2. Subsection 4.3 outlines viability approach and describes a grid numerical method for constructing leadership kernels, i.e. maximal domains where the control can keep the dynamic system infinitely long for all unpredictable admissible disturbances. The idea of the proof of robustness of the heuristic controller is then explained in terms of leadership sets. In Subsection 4.4, a modified grid algorithm for the treatment of not quickly changing disturbances is given.
Section 5 describes results of numerical simulations. It is shown that the controller introduced in Section 4 is able to correct the impaired autoregulation, i.e. the controller is robust against all unpredictable admissible not quick changes of MAP, pCO2, and CVP.
In this Section, a model of cerebral autoregulation for preterm infants, originally introduced in [14], is recalled and extended to account for effects of carbon dioxide vasoreactivity.
A hierarchical model of cerebrovascular network proposed in [17] is used to compute CBF. It is supposed that 9 arterial, one capillary, and 9 venous compartments are sequentially connected. In the case of Poiseuille flow of a Newtonian liquid, CBF can be evaluated using Kirchhoff's law as follows (see [21] for details):
CBF=(pa−pv)(19∑i=18μℓiπmir4i)−1. | (2.1) |
Here, pa=MAP is the mean arterial pressure, pv=CVP is the cerebral venous pressure, ℓi and ri are the length and radius of vessels of the ith compartment, respectively, and μ is the dynamic viscosity of blood.
To describe the process of autoregulation, assume that the vascular radius ri depends on the arterial pressure pa as follows:
ri=r∗i[(p∗a−pv)/(pa−pv)]1/4, |
where p∗a is a baseline value of pa, and r∗i is a reference value of ri corresponding to p∗a. It is easy to see that such a choice stabilizes CBF because ri appears as the fourth power in formula (2.1) for CBF.
To fit this model to experimental data on cerebral autoregulation in preterm infants (see [3,22,23,24], the following modification is introduced (see also [14]):
ri=r∗i[(p∗a−pv)/(pa−pv)]1/4P(pa,a1,a2,a3), | (2.2) |
P(pa,a1,a2,a3)=1+a1(pa−p∗a)+a2(pa−p∗a)2+a3(pa−p∗a)3. | (2.3) |
Therefore,
CBF(pa,pv,a1,a2,a3)=(pa−pv)kage(19∑i=18μℓiπmir4i)−1, | (2.4) |
where ri are given by formulas (2.2) and (2.3), and kage is a scale factor to adjust CBF to the age of infants. The values p∗a=35 [mmHg] and kage=0.08 correspond to hemodynamic system of premature infants of 31-34 weeks' gestational age, the brain weight of 260 g and the CBF of 15.5 mL/100 g/min (cf. [22,24]). The value of pv is set to be equal to 5 [mmHg] in the fitting procedure. The values of ℓi, mi, and r∗i, i=1,...,19, are taken from [17].
The coefficients a1,a2, and a3 are fitted through the minimization of the residual
R(a1,a2,a3)=5∑k=1[CBF(pka,pv,a1,a2,a3)−CBFk]2, | (2.5) |
where the pairs (pka[mmHg],CBFk[mL/s]),k=1,...5, are chosen according to experimental data from [22,24] as follows:
(κ⋅20,0.24);(κ⋅30,0.67);(κ⋅34,0.67);(κ⋅38,0.67);(κ⋅50,2). |
Here κ=133.322 Pa/mmHg is the conversion factor from mmHg to Pa.
Note that the polynomial (2.3) of degree 3 and the five data points are sufficient to approximate the autoregulation plateau experimentally found in [22].
The values of the minimizers of (2.5) read
a1=8.384e-7 Pa−1,a2=5.718e-9 Pa−2,a3=3.518e-11 Pa−3. |
Although the coefficients ai are very small, the last three summands in (2.3) are not too small because the pressure is measured in pascals.
To account for dilation/constriction of blood vessels caused by the change of carbon dioxide partial pressure in arterial blood, the following dependence of the reference radii r∗i,i=1,...,19, on the partial CO2 pressure, pCO2, is adopted from [17]:
r∗i=r∗i(pCO2)=r0i⋅(1+ci⋅(pCO2−pCO02)), | (2.6) |
where ci is the pCO2 reactivity coefficient, pCO02 a baseline value of the partial CO2 pressure, and r0i the baseline vessel radius.
In this section, the autoregulation model from Sect. 2 is modified to describe autoregulation defects in preterm infants. Namely, it will be supposed that the horizontal plateau observed in Figure 1 may be distorted due to the change in pCO2 and CVP, which also extends the model of impaired cerebral autoregulation introduced by the authors in [16]. Making the substitutions
pa←x1,pCO2←x2,a1←ϕ(x2)a1,r∗i←(1+x3)r∗i,pv←x4, |
in formulas (2.2)-(2.4) and (2.6) yields a function
q(x1,x2,x3,x4)=kage(x1−x4)(∑19i=18μℓiπmir4i(x1,x2,x3,x4))−1,ri(x1,x2,x3,x4)=(1+x3)r0i(1+ci(x2−pCO02))×[(p∗a−x4)/(x1−x4)]1/4(ϕ(x2)a1(x1−p∗a)+a2(x1−p∗a)2+a3(x1−p∗a)3) |
describing the disturbed autoregulation. Here; x1 denotes the mean arterial pressure, pa; x2 denotes the partial CO2 pressure, pCO2; the multiplier ϕ(x2) changes the slope of the plateau observed in Figure 1; the variable x3 modifies the vascular volume to restore the horizontal autoregulation plateau; and x4 stands for the venous pressure.
The dependence ϕ=ϕ(x2) is reconstructed from the experimental data for preterm infants (see [3] where the slope s(pCO2) of the curve presenting CBF velocity versus pCO2 is measured). The function ϕ is found numerically from the relation
CBF(¯pa,p∗v,ϕ(x2)a1,a2,a3)−CBF(pa_,p∗v,ϕ(x2)a1,a2,a3)¯pa−pa_=S⋅s(x2), |
where pa_ and ¯pa are the starting and finishing values of the arterial blood pressure for the autoregulation plateau, p∗v=5 mmHg is the same value used in the fitting procedure for the coefficients a1, a2, a3, and S is the estimate of the mean cross-section area of the blood vessels. The graph of ϕ is shown in Figure 2. Obviously, the relation ϕ(pCO02)=1 holds.
Note that the function q(x1,pCO02,0,p∗v) represents the normal autoregulation response shown in Figure 1. Therefore, the function
ω(x1,x2,x3,x4):=q(x1,x2,x3,x4)−q(x1,pCO02,0,p∗v) | (3.1) |
can be considered as the measure of the violation of autoregulation. This discrepancy will be used below to define a state constraint in the forthcoming conflict control system.
In this section, a conflict control model governing the variables x1,x2,x3,x4 and accounting for different state constraints will be stated. Moreover, a feedback control based on the discrepancy between the current and nominal CBF will be designed.
Consider now the following conflict control problem (differential game [20]):
˙x1=−k1(x1−v1),˙x2=−k2/ϕ′(x2)(ϕ(x2)−v2),˙x3=−k3(x3−u),˙x4=−k4(x4−v3) | (4.1) |
where u is the control variable of the first player, which can be interpreted as input of a medicament that dilates or constricts blood vessels with some time delay defined by the coefficient k3. The disturbance variables v1, v2, and v3 are at the disposal of the second player. They have effect on the arterial pressure, x1, the partial CO2 pressure, x2, and the venous pressure, x4, respectively. The coefficients k1, k3, and k4 define the corresponding time delays. It should be noted that the second equation is equivalent to the following: ˙χ2=−k2(χ2−v2) with χ2=ϕ(x2). Thus, the model (4.1) uses the physical variable x2:=pCO2 instead of the abstract parameter χ2 utilized (up to notation) in the model presented in [16].
The state variables are assumed to be constrained as follows:
x1∈[28,42]mmHg,x2∈[35,60]mmHg,x3∈[−0.2,0.2],x4∈[0,10]mmHg. | (4.2) |
The disturbances and control are constrained as follows:
v1∈[28,42]mmHg,v2∈[30,55]mmHg,u∈[−0.2,0.2],v3∈[0,10]mmHg. | (4.3) |
The right-hand side values of the constraints in (4.2) and (4.3) are chosen based on clinical measurements presented in [23,25,26]. The bounds on control inputs are set to enable the control to compensate large deviations of CBF from the autoregulation plateau within 10 minutes.
It should be stressed that equations (4.1) represent the so called PT1 filters, where the coefficients k1, k2, k3, and k4 define the time constants, reaction times to input changes. It is assumed that time in system (4.1) is measured in minutes, and the coefficients are chosen as follows:
k1=0.01min−1, k2=0.01min−1, k3=0.1min−1, k4=0.01min−1. |
So, we suppose that the reaction time to changes of the arterial pressure, partial carbon dioxide pressure, and cerebral venous pressure is about 100 minutes, which is consistent with our clinical observations. Moreover, the response to the medication occurs in about 10 minutes after injection.
The dynamical system (4.1)–(4.3) will be considered as an antagonistic differential game with two opposite players. The objective of the first player is to ensure the constraints (4.2) and the following constraint:
|ω(x1(t),x2(t),x3(t),x4(t))|≤ω0,t∈[0,∞), | (4.4) |
where ω0=2mL/min is the admissible deviation from the perfect autoregulation performance.
Note that the highly nonlinear state constraint (4.4) is the central part of the model, since it reflects the deviation of computed CBF from the reference value, whereas equations (4.1) only restrict the rate of change of model variables in time. Obviously, the constraints (4.2) and (4.4) should hold for all disturbances v1(⋅), v2(⋅), and v3(⋅) satisfying the bounds (4.3). It is easy to prove that the state constraints (4.2) hold whenever the bounds (4.3) hold. Thus, the control u should not spend much trouble on keeping these state constraints. The second player strives to violate the constraint (4.4).
It is easy to check that the following relation holds:
q(x1,x2,x3,x4)=(1+x3)4⋅q(x1,x2,0,x4) | (4.5) |
According to the relation (4.5) and the definition (3.1), it is sufficient to set
x3=(q(x1,pCO02,0,p∗v)q(x1,x2,0,x4))14−1 | (4.6) |
to provide ω=0.
Assume now that the control u is chosen at a time instant t and at a state {x1=x1(t), x2=x2(t), x3=x3(t), x4=x4(t)} as follows:
u=(q(x1,pCO02,0,p∗v)q(x1,x2,0,x4))14−1. | (4.7) |
Due to the third equation of (4.1), it should be expected that the variable x3(t) will follow the control u(t) because the coefficient k3 is essentially larger than k1,k2, and k4. Thus, the relation (4.6) will be approximately fulfilled, which will provide the condition ω≈0.
Taking into account the relation (4.5) and using the notation x=(x1,x2,x3,x4), we can rewrite (4.7) as follows:
u(x)=(q(x1,pCO02,0,p∗v)q(x1,x2,x3,x4))14(1+x3)−1. | (4.8) |
Here, q(x1,x2,x3,x4) represents the current CBF that can physically be measured, and q(x1,pCO02,0,p∗v) is the reference CBF at arterial pressure x1, which is supposed to be known. The variable x3 depends on the amount of medicine injected, and, therefore, the current value of x3 can be easily computed.
Note that the control (4.8) compensates the deviation of CBF from its reference value not perfectly because of a time delay caused by the dynamics of the third equation of (4.1). This will be seen below from simulations with a discrete control scheme. To improve the quality of the control (4.8), it is reasonable to permit an oversteering that is proportional to the discrepancy ω(x1,x2,x3,x4) with the experimentally obtained coefficient of 0.9. Thus, we arrive at a modified control
u∗(x)=(q(x1,pCO02,0,p∗v)−0.9⋅ω(x1,x2,x3,x4)q(x1,x2,x3,x4))14(1+x3)−1. |
The role of the oversteering will be clearly seen in simulations based on discrete-time control scheme with large time sampling intervals. Finally, the following feedback control is suggested:
˜u(x)={−μ,u∗(x)<−μ,μ,u∗(x)>μ,u∗,u∗(x)∈[−μ,μ], μ=0.2. | (4.9) |
Thus, to compute the control (intake of medicine), the current CBF and the arterial blood pressure are measured, the reference CBF value for the actual arterial pressure is calculated by the formula, and the accumulated value of injected medicine is evaluated. We will see that the control (4.9) is able to better track the reference CBF curve, i.e., provide smaller values of the discrepancy ω. Moreover, using a relaxed algorithm for computing leadership sets, we will prove that ˜u guarantees keeping trajectories within the state constraint (4.2) for all admissible disturbances with not quick variations. For consistency, a base algorithm for constructing leadership kernels will be sketched in the next section.
The notion of viability kernel (see [19]) is used for control problems, whereas, in the case of differential games, the notions of leadership and discriminating kernels are conventionally utilized.
Leadership kernel is the maximal subset of the state constraint where the system trajectories can remain arbitrary long if the first player utilizes an appropriate pure feedback control, and the second player generates any admissible disturbances (see [27,28]). In contrast, the notion of discriminating kernel corresponds to the case where the first player can exactly measure current actions of the second one to use the so-called counter feedback strategies (see[20]). Obviously, this case is not realistic.
In this section, viability approach and a numerical method developed in [21,29] for constructing leadership kernels will be briefly outlined. Note that in our case leadership and discriminating kernels coincide, because Isaacs' saddle point condition holds for dynamic system (4.1).
Remember that x=(x1,x2,x3,x4) is the state vector and denote by f(x,u,v), v=(v1,v2,v3), the right-hand side of system (4.1). Let P and Q be compact sets defining the bounds (4.3) on the control and disturbances, respectively.
The results of this subsection hold for arbitrary dimension n of the state space and rather general function f (see e.g. [29]), which is reflected in the description below.
Let u→v(u) be a Borel measurable function with values in Q. Consider the differential inclusion
˙x∈Fv(⋅)(x)=¯co{f:f(x,u,v(u)),u∈P}. | (4.10) |
Definition 4.1 (leadership property). A set K⊂Rn has leadership property if for any x∗∈K and any Borel measurable function v(⋅) with values in Q there exists a solution x(⋅) of (4.10) such that x(0)=x∗ and x(t)∈K for all t≥0.
Definition 4.2 (leadership kernel). For a given compact set G⊂Rn denote by Lead(G) the largest subset of G with the leadership property. This subset is called the leadership kernel of G.
Let
Gλ={x∈Rn, g(x)≤λ} |
be a family of state constraints. In the simulations below, g(x)=maxigi(x) is set, where g1(x)=|x1−35|−7, g2(x)=|x2−47.5|−12.5, g3(x)=|x3|−0.2, g4(x)=|x4−5|−5, and g5(x)=|ω(x)|−ω0. The first four functions account for the state constraints (4.2) and the last one is responsible for the constraint (4.4). Therefore, Gλ is equivalent to the constraints (4.2) and (4.4) if λ=0.
It is required to construct a function V, such that
Lead(Gλ)={x∈Rn, V(x)≤λ}. | (4.11) |
Denote by δ>0 the time step length. Let h:=(h1,...,hn) be the space grid steps with |h|=max1≤i≤nhi. For any continuous function V:Rn→R introduce the following upwind operator:
Π(V;δ,h)(x)=V(x)+δminu∈Pmaxv∈Qn∑i=1(prightif+i+pleftif−i), |
where fi are the components of f, and
a+=max{a,0},a−=min{a,0},prighti=[V(x1,...,xi+hi,...,xn)−V(x1,...,xi,...,xn)]/hi,plefti=[V(x1,...,xi,...,xn)−V(x1,...,xi−hi,...,xn)]/hi. |
Note that prighti and plefti are typical finite differences in Godunov-like methods for hyperbolic problems. The operator Π will be applied to grid functions to return grid functions.
Denote Vh(xi1,...,xin)=V(i1h1,...,inhn), gh(xi1,...,xin)=g(i1h1,...,inhn) being the restrictions of V and g to the grid. Let a sequence {δℓ} is chosen, such that δℓ→0, ∑∞ℓ=0δℓ=∞.
Consider the following grid scheme:
Vhℓ+1=max{Π(Vhℓ;δ,h),gh}, Vh0=gh,ℓ=0,1,…. | (4.12) |
It can be proven that Vhℓ monotonically point-wise converges to a grid function Vδ,h, and this function approximates the function V defining leadership kernels according to formula (4.11) if δ and h are sufficiently small.
Proposition 1 (see [30] for the proof). Let |Fv(x)|≤B for all x∈Rn and δℓ/|h|≤1/(B√n) for all ℓ. Then Vhℓ↗Vh as ℓ→∞.
Note that in the case n=4 the relation δ/|h|≤(2B)−1 should hold. Other secondary stability conditions can be found in [31,21].
The estimate |Vh−V|≤C√|h| is expected (cf. [31,21]), and therefore Vhℓ approximates V if ℓ is large and |h| is small. The stopping criterion in the computation is sup grid|Vhℓ+1−Vhℓ|≤ε.
To be sure that the control (4.9) works well for all unpredictable admissible disturbances that change not quickly, the following technique can be applied. First, extend the system (4.1) by adding three PT1 filters for v1,v2, and v3 with the same time constant k=0.1 min−1, which will restrict the rate of change of the disturbances. Thus, the extension of (4.1) looks as follows:
˙vi=−k(vi−ˉvi), i=1,2,3, | (4.13) |
where ˉvi,i=1,2,3 , are new disturbances constrained in the same way as the old ones, cf. (4.3).
Second, put ˜u into the extended system (4.1) and compute a leadership set of the state constraints (4.2), (4.4) according to a method described below. If this leadership set is nonempty, the heuristic rule (4.9) is appropriate to work against all unpredictable admissible disturbances that change not quickly. Note, that the leadership set will be computed in the four-dimensional space of the variables x1,x2,x3,x4. The variables v1,v2, and v3 are not considered as new states, but as accumulated optimal values of the new artificial controls.
Let wi,i=1,2,3, be given grid functions. Consider now the following operator:
˜Π(V;δ,h,w1,w2,w3)(x)=V(x)+δmax(ˉv1,ˉv2,ˉv3)4∑i=1(prighti˜f+i+plefti˜f−i), | (4.14) |
where
˜f1(x,ˉv1)=−k1(x1−v1),˜f2(x,ˉv2)=−k2/ϕ′(x2)(ϕ(x2)−v2),˜f3(x)=−k3(x3−u),˜f4(x,ˉv3)=−k4(x4−v3) |
with the substitution
v1=w1(x)−kδ(w1(x)−ˉv1), | (4.15) |
v2=w2(x)−kδ(w2(x)−ˉv2), | (4.16) |
v3=w3(x)−kδ(w3(x)−ˉv3),u=˜u(x). | (4.17) |
The recurrent computation scheme (4.12) transforms now into the following one:
Vh0(x)=gh(x),w01(x)=w02(x)=w03(x)=0, | (4.18) |
Vhℓ+1(x)=max{˜Π(Vhℓ;δ,h,wℓ1,wℓ2,wℓ3)(x),gh(x)}, | (4.19) |
wℓ+1i(x)=wℓi(x)+kδ(wℓi(x)−ˉvℓ,maxi(x)),i=1,2,3. | (4.20) |
Here, ˉvℓ,maxi(x),i=1,2,3 , are the maximizes, over ˉvi,i=1,2,3 , when computing the term
˜Π(Vhℓ;δ,h,wℓ1,wℓ2,wℓ3)(x) |
according to formula (4.14).
Note that the relations (4.15)–(4.17) represent, at each node x, a time-step approximation of equations (4.13) with initial conditions wi(x),i=1,2,3. These initial conditions are obtained from the previous time step as solutions of (4.13) with some optimal values of ˉvi,i=1,2,3. Therefore, wℓi(x),i=1,2,3, accumulate, according to equation (4.13), values of ˉvi,i=1,2,3, that are optimal at the grid point x, see (4.20). Taking into account that first player's control is a prescribed function, the relation (4.20) expresses the dynamic programming optimality principle. It should be noted that the starting functions w0i(x),i=1,2,3, can be chosen arbitrary, under the constraints (4.3) with vi being replaced by w0i(x) for i=1,2,3.
It can be proven that the modified grid scheme (4.18)–(4.20) is monotone and a similar convergence result as in Proposition 1 can be established. Moreover, the limiting function does not increase on the trajectories of the system (4.1), (4.13), with u=˜u(x), for all admissible controls ˉvi(t),i=1,2,3. Therefore, level sets of the limiting function are leadership sets.
Here, simulation results for system (4.1) will be presented, and the robustness of the above proposed heuristic feedback controller (4.9) against not quickly changing admissible disturbances v1,v2,v3 will be checked.
First, the ability of the control ˜u against open-loop disturbances, such as e.g. those shown in Figure 3, is tested. Here, the behavior of the mean arterial pressure is as follows: It increases from 28 mmHg to 42 mmHg during first 16 hours, remains constant (equals 42 mmHg) during further 17 hours, linearly decreases to 28 mmHg during 1.5 days, then remains constant (equals 28 mmHg) during 2 days, and finally increases to 42 mmHg during 10 hours. The pCO2 value changes periodically every 30 hours, taking values between 30 and 55 mmHg. For the cerebral venous pressure, a sinusoidal change with the amplitude of 5 mmHg and the periodicity of approximately 31 hours is set. The control (4.9) is applied continuously. Note that the chosen disturbances obey the constraints (4.3) and simulate apparent clinical conditions.
In Figure 4, the resulting autoregulation curve (red) computed with the control (4.9) and above described disturbances is shown. Additionally, the reference autoregulation response curve (green) is plotted. Note that the graphs are plotted parametrically, i.e., x=x1(t), y=q(x1(t),x2(t),x3(t),x4(t)) and x=x1(t), y=q(x1(t),pCO02,0,p∗v), respectively. The time development of the feedback control (4.9) is presented in Figure 5 in green color. In red color, the time development of the variable x3 is shown. One can see a good closeness of x3(t) and u(t) in the case of continuous-time control scheme.
A discrete-time control scheme with 2 hours time sampling is tested. This means that the feedback control is computed every 2 hours according to the rule (4.9) and is kept constant in between. It can be interpreted as a regular intake of a medicine like e.g. indomethacin [32] for lowering CBF. The resulting autoregulation curve for such a discrete feedback control and disturbances shown in Figure 3 is presented in Figure 6. One can see that, even in this less favorable case, deflections of CBF do not exceed 2 mL/min, which is visualized with the help of a strip around the undisturbed autoregulation curve. The time development of the discrete feedback control with oversteering (green line) is depicted in Figure 7. Simultaneously the time development of the variable x3 (red line) is shown there. It is seen that these lines remain close to each other, as it was in the case of continuous-time control scheme. It should be noted that the control without oversteering term is not able to sufficiently push the trajectory to the reference curve because of a time delay caused by the dynamics of the third equation of (4.1), which produces an unsatisfying result (see Figure 8).
Note that the oscillations of the CBF curve around the autoregulation plateau, observable in Figures 6 and 8, are caused by the discrete-time control scheme with large sampling interval of 2 hours. During the sampling time, the disturbance pushes the CBF curve away from the autoregulation plateau, and the controller should bring the trajectory back to the plateau level. Apparently, it might be possible to modify the control in order to smooth the cusps, which, however, will not decrease the amplitude of the oscillations. Moreover, these oscillations are not dangerous because they lie in admissible physiological range (5% of the CBF plateau value).
To prove that the proposed feedback control is powerful against all "slow" admissible disturbances constrained similar to (4.3), this control is inserted into the extended system (4.1) and (4.13), and the leadership kernel of the state constraints (4.2), (4.4) is computed as described in Subsection 4.4. Note that the nonlinearity both of the system (4.1) and of the state constraint (4.4) as well as the relatively high dimension of the state vector (four) make the computation of the leadership kernel to be a complex problem.
The computation was performed on the SuperMUC system at the Leibniz Supercomputing Centre of the Bavarian Academy of Sciences and Humanities. The problem was parallelized between 25 compute nodes with 16 cores per node. A grid of 100×100×100×100 cells was used, and about 104 steps of the algorithm (4.18)–(4.20) were performed. The runtime was about 1 hour.
The result of the computation is presented in Figure 9. The section of the leadership set at x4=p∗v is shown. The vertical axis measures the discrepancy ω. The non-emptiness of the computed leadership kernel means that the control (4.9) is able to keep trajectories within the state constraints (4.2) and (4.4) for all admissible disturbances with not quick variations.
The scientific awareness about possible controlling of impaired cerebral autoregulation in preterm infants is still far from satisfactory. To the best of our knowledge, there are no mathematical models available to control impaired cerebral autoregulation in preterm infants. Adequate mathematical modeling can help to improve understanding of cerebral circulation in premature brain and assess the influence of the most important factors like MAP, pCO2, and CVP on the cerebral blood flow.
In the current paper, impaired cerebral autoregulation in preterm infants has been mathematically modeled. Feedback control that can be interpreted as a medication strategy preventing large deviations of CBF from some physiological (reference) autoregulation curve has been developed. The slope of the autoregulation plateau due to hypercapnia is reconstructed from the experimental data for preterm infants in [3]. The ranges of change for mean arterial pressure, partial CO2 pressure, and cerebral venous pressure are consistent with clinical measurements in [26], however, because of the lack of measurements, the time-evolutions of these variables used in our simulations are heuristic curves that satisfy the chosen constraints.
The model developed can be a reasonable starting point to improve understanding of autoregulation in preterm infants. Apparently, such a modeling may support refining clinical therapies and monitoring strategies. Viability theory is applied to prove that the proposed control works well against all admissible disturbances of arterial, venous, and partial CO2 pressure. Moreover, satisfactory result is obtained in the case of discrete-time control scheme with time sampling of 2 hours. The advantage of the controller is its simple usage based on only two measurements, CBF and MAP, at each sampling time instant. However, it would be interesting to improve the controller, for example, to avoid oversteering, by using nonlinear automatic control techniques such as output regulation or disturbance rejection.
It is worth to mention that the described viability approach and numerical method can be used in various applications to assess the quality of controls. Future work will be focused on coupling of the control-based models of cerebral autoregulation with models of cerebral vascular systems accounting for peculiarities of the immature brain (see e.g. [33]).
The authors gratefully acknowledge financial support of the Klaus Tschira Foundation, W{ü}rth Foundation, and Buhl-Strohmaier-Foundation. The last author appreciates the Markus Würth Professorship at the Technical University of Munich. Computer resources have been provided by the Gauss Centre for Supercomputing/Leibniz Supercomputing Centre under grant: pr74lu.
The authors declare that there is no conflict of interest regarding the publication of this article.
[1] | J. Donnelly, M. J. Aries and M. Czosnyka, Further understanding of cerebral autoregulation at the bedside: possible implications for future therapy, Expert Rev. Neurother., 15 (2015), 169–185. |
[2] | M. van de Bor and F. J. Walther, Cerebral blood flow velocity regulation in preterm infants, Biol. Neonate, 59 (1991), 329–335. |
[3] | J. R. Kaiser, C. H. Gauss and D. K.Williams, The effects of hypercapnia on cerebral autoregulation in ventilated very low birth weight infants, Pediatr. Res., 58 (2005), 931–935. |
[4] | J. J. Volpe, Intracranial hemorrhage: germinal matrix hemorrhage of the premature infant, in Neurology of the Newborn, 48 (2001), WB Saunders, Philadelphia, pp. 435. |
[5] | T. Lekic, D. Klebe, R. Poblete, et al., Neonatal brain hemorrhage (NBH) of prematurity: translational mechanisms of the vascular-neural network, Curr. Med. Chem., 22 (2015), 1214– 1238. |
[6] | O. Pryds and A. D. Edwards, Cerebral blood flow in the newborn infant, Arch. Dis. Child- Fetal, 74 (1996), F63–F69. |
[7] | P. Ballabh, Intraventricular hemorrhage in premature infants: mechanism of disease, Pediatr. Res., 67 (2010), 1–8. |
[8] | M. Ursino and C. A. Lodi, A simple mathematical model of the interaction between intracranial pressure and cerebral hemodynamics, J. Appl. Physiol., 82 (1997), 1256–1269. |
[9] | R. B. Panerai, S. L. Dawson and J. F. Potter, Linear and nonlinear analysis of human dynamic cerebral autoregulation, Am. J. Physiol.- Heart C, 277 (1999), H1089–H1099. |
[10] | M. S. Olufsen, A. Nadim and L. A. Lipsitz, Dynamics of cerebral blood flow regulation explained using a lumped parameter model, Am. J. Physiol.- Reg. I, 282 (2002), R611–R622. |
[11] | J. Alastruey, S. M. Moore, K.H. Parker, et al., Reduced modeling of blood flow in cerebral circulation: Coupling 1-D, 0-D and cerebral auto-regulation models, Int. J. Numer. Meth. Fl., 56 (2008), 1061–1067. |
[12] | V. Z. Marmarelis, D. C. Shin and R. Zhang, Linear and nonlinear modeling of cerebral flow autoregulation using principal dynamic modes, Open Biomed. Eng. J., 6 (2012), 42–55. |
[13] | B. Spronck, E. G. Martens, E. D. Gommer, et al., A lumped parameter model of cerebral blood flow control combining cerebral autoregulation and neurovascular coupling, Am. J. Physiol.- Heart C, 303 (2012), H1143–11H53. |
[14] | R. Lampe, N. Botkin, V. Turova, et al., Mathematical modelling of cerebral blood circulation and cerebral autoregulation: towards preventing intracranial hemorrhages in preterm newborns, Comput. Math. Method M., (2014), 1–9. |
[15] | S. Payne, Cerebral Autoregulation: Control of Blood Flow in the Brain, Springer, Berlin, 2016. |
[16] | N. Botkin, V. Turova and R. Lampe, Feedback control of impaired cerebral autoregulation in preterm infants: Mathematical modelling, in 25th Mediterranean Conference on Control and Automation (MED), IEEE, (2017), 229–234. |
[17] | S. K. Piechnik, P. A. Chiarelli and P. Jezzard, Modelling vascular reactivity to investigate the basis of the relationship between cerebral blood volume and flow under CO2 manipulation, NeuroImage, 39 (2008), 107–118. |
[18] | U. H. Thome and N. Ambalavanan, Permissive hypercapnia to decrease lung injury in ventilated preterm neonates, Semin. Fetal Neonat. M., 14 (2009), 21–27. |
[19] | J. P. Aubin, Viability Theory, Birkhäuser, Boston, 1991. |
[20] | N. N. Krasovskii and A. I. Subbotin, Game-theoretical control problems, Springer, New York, 1988. |
[21] | N. D. Botkin and V. L. Turova, Numerical construction of viable sets for autonomous conflict control systems, Mathematics, 2 (2014), 68–82. |
[22] | J. L. LeFlore and W. D. Engle, Clinical factors influencing blood pressure in the neonate, NeoReviews, 3 (2002), e145–e150. |
[23] | G. D. T. Inglis, K. R. Dunster and M. W. Davies, Establishing normal values of central venous pressure in very low birth weight infants, Physiol. Meas., 28 (2007), 1283. |
[24] | J. M. Perlman, Neurology: Neonatology questions and controversies, Saunders/Elsevier, Philadelphia, 2008. |
[25] | T.A. Polovova, Hemodynamics Features with Different Methods of Respiratory Support in Children with Extremely Low Body Weight, Ph.D thesis, Ural Scientific Research Institute of Maternity and Infancy Protection of the Ministry of Health of the Russian Federation in Ekaterinburg, 2014 (in Russian). |
[26] | R. Lampe, V. Turova, N. Botkin, et al., Postnatal paraclinical parameters associated to occurrence of intracerebral hemorrhage in preterm infants, Neuropediatrics, (2018). |
[27] | P. Cardaliaguet, A differential game with two players and one target, SIAM J. Control Optim., 34 (1996), 1441–1460. |
[28] | P. Cardaliaguet, M. Quincampoix and P. Saint-Pierre, Set valued numerical analysis for optimal control and differential games, in M. Bardi, T. E. S. Raghavan, T. Parthasarathy (eds.) Stochastic and Differential Games: Theory and Numerical Methods. Annals of the International Society of Dynamic Games 4 (1999), Birkhäuser, Boston, 177–274. |
[29] | N. Botkin, V. Turova, J. Diepolder, et al., Aircraft control during cruise flight in windshear conditions: viability approach, Dyn. Games Appl., (2017), 1–15. |
[30] | N. D. Botkin and V. L. Turova, Examples of computed viability kernels, Trudy Inst. Mat. i Mekh. UrO RAN, 21 (2015), 306–319. |
[31] | N. D. Botkin, K.-H. Hoffmann, N. Mayer, et al., Approximation schemes for solving disturbed control problems with non-terminal time and state constraints, Analysis, 31 (2011), 355–379. |
[32] | A. D. Edwards, J. S. Wyatt, C. Richardson, et al., Effects of indomethacin on cerebral haemodynamics in very preterm infants, The Lancet, 335 (1990), 1491–1495. |
[33] | N. D. Botkin, A. E. Kovtanyuk, V. L. Turova, et al., Direct modeling of blood flow through the vascular network of the germinal matrix, Comput. Biol. Med., 92 (2018), 147–155. |
1. | Elisabeth M.W. Kooi, Anne E. Richter, Cerebral Autoregulation in Sick Infants, 2020, 47, 00955108, 449, 10.1016/j.clp.2020.05.003 | |
2. | Stephen A. Back, Joseph J. Volpe, 2025, 9780443105135, 523, 10.1016/B978-0-443-10513-5.00019-X |