Linking state-and-transition simulation and timber supply models for forest biomass production scenarios

1.   Introduction

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 $ CO_2 $ level ($ pCO_2 $) ^[3]. Increased $ pCO_2 $ 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 $, $ pCO_2 $, 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 $ pCO_2 $ 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 $, $ pCO_2 $, 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 $, $ pCO_2 $, 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 $, $ pCO_2 $, and $ CVP $.

2.   Mathematical model of cerebral autoregulation

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):

Here, $ p_a = MAP $ is the mean arterial pressure, $ p_v = CVP $ is the cerebral venous pressure, $ \ell_i $ and $ r_i $ are the length and radius of vessels of the $ i $th compartment, respectively, and $ \mu $ is the dynamic viscosity of blood.

To describe the process of autoregulation, assume that the vascular radius $ r_i $ depends on the arterial pressure $ p_a $ as follows:

where $ p_a^* $ is a baseline value of $ p_a $, and $ r^*_i $ is a reference value of $ r_i $ corresponding to $ p_a^* $. It is easy to see that such a choice stabilizes $ CBF $ because $ r_i $ 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]):

Therefore,

where $ r_i $ are given by formulas (2.2) and (2.3), and $ k_{age} $ is a scale factor to adjust $ CBF $ to the age of infants. The values $ p^*_a = 35$ [mmHg] and $ k_{age} = 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 $ p_v $ is set to be equal to 5 [mmHg] in the fitting procedure. The values of $ \ell_i $, $ m_i $, and $ r^*_i $, $ i = 1, ..., 19 $, are taken from ^[17].

The coefficients $ a_1, a_2 $, and $ a_3 $ are fitted through the minimization of the residual

where the pairs $ \big(p^k_a\; {\rm[mmHg]}, CBF^k\; {\rm[mL/s]}), k = 1, ...5, $ are chosen according to experimental data from ^[22,24] as follows:

Here $ \kappa = 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

Although the coefficients $ a_i $ 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 $ CO_2 $ pressure, $ pCO_2 $, is adopted from ^[17]:

where $ c_i $ is the $ pCO_2 $ reactivity coefficient, $ pCO_2^0 $ a baseline value of the partial $ CO_2 $ pressure, and $ r_i^0 $ the baseline vessel radius.

3.   Modeling impaired cerebral autoregulation

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 $ pCO_2 $ and $ CVP $, which also extends the model of impaired cerebral autoregulation introduced by the authors in ^[16]. Making the substitutions

in formulas (2.2)-(2.4) and (2.6) yields a function

describing the disturbed autoregulation. Here; $ x_1 $ denotes the mean arterial pressure, $ p_a $; $ x_2 $ denotes the partial $ CO_2 $ pressure, $ pCO_2 $; the multiplier $ \phi(x_2) $ changes the slope of the plateau observed in Figure 1; the variable $ x_3 $ modifies the vascular volume to restore the horizontal autoregulation plateau; and $ x_4 $ stands for the venous pressure.

The dependence $ {\phi} = \phi(x_2) $ is reconstructed from the experimental data for preterm infants (see ^[3] where the slope $ s(pCO_2) $ of the curve presenting $ CBF $ velocity versus $ pCO_2 $ is measured). The function $ \phi $ is found numerically from the relation

where $ \underline{p_{a}} $ and $ \overline{p_a} $ 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 $ a_1 $, $ a_2 $, $ a_3 $, and $ {\mathcal S} $ is the estimate of the mean cross-section area of the blood vessels. The graph of $ \phi $ is shown in Figure 2. Obviously, the relation $ \phi(pCO_2^0) = 1 $ holds.

Note that the function $ q(x_1, pCO_2^0, 0, p_v^*) $ represents the normal autoregulation response shown in Figure 1. Therefore, the function

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.

4.   Conflict control problem and heuristic quasi-optimal control

In this section, a conflict control model governing the variables $ x_1, x_2, x_3, x_4 $ 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.

4.1.   Conflict control problem

Consider now the following conflict control problem (differential game ^[20]):

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 $ k_3 $. The disturbance variables $ v_1 $, $ v_2 $, and $ v_3 $ are at the disposal of the second player. They have effect on the arterial pressure, $ x_1 $, the partial $ CO_2 $ pressure, $ x_2 $, and the venous pressure, $ x_4 $, respectively. The coefficients $ k_1 $, $ k_3 $, and $ k_4 $ define the corresponding time delays. It should be noted that the second equation is equivalent to the following: $ \dot{\chi}_2 = -k_2 (\chi_2-v_2) $ with $ \chi_2 = \phi(x_2) $. Thus, the model (4.1) uses the physical variable $ x_2: = pCO_2 $ instead of the abstract parameter $ \chi_2 $ utilized (up to notation) in the model presented in ^[16].

The state variables are assumed to be constrained as follows:

The disturbances and control are constrained as follows:

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 $ k_1 $, $ k_2 $, $ k_3 $, and $ k_4 $ 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:

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:

where $ \omega_0 = 2\; \text{mL/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 $ v_1(\cdot) $, $ v_2(\cdot) $, and $ v_3(\cdot) $ 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).

4.2.   Control design

It is easy to check that the following relation holds:

According to the relation (4.5) and the definition (3.1), it is sufficient to set

to provide $ \omega = 0 $.

Assume now that the control $ u $ is chosen at a time instant $ t $ and at a state $ \big\{x_1 = x_1(t), \ x_2 = x_2(t), \ x_3 = x_3(t), \ x_4 = x_4(t)\big\} $ as follows:

Due to the third equation of (4.1), it should be expected that the variable $ x_3(t) $ will follow the control $ u(t) $ because the coefficient $ k_3 $ is essentially larger than $ k_1, k_2 $, and $ k_4 $. Thus, the relation (4.6) will be approximately fulfilled, which will provide the condition $ \omega\approx 0 $.

Taking into account the relation (4.5) and using the notation $ x = (x_1, x_2, x_3, x_4) $, we can rewrite (4.7) as follows:

Here, $ q(x_1, x_2, x_3, x_4) $ represents the current $ CBF $ that can physically be measured, and $ q(x_1, pCO_2^0, 0, p_v^*) $ is the reference $ CBF $ at arterial pressure $ x_1 $, which is supposed to be known. The variable $ x_3 $ depends on the amount of medicine injected, and, therefore, the current value of $ x_3 $ 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 $ \omega(x_1, x_2, x_3, x_4) $ with the experimentally obtained coefficient of 0.9. Thus, we arrive at a modified control

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:

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 $ \omega $. Moreover, using a relaxed algorithm for computing leadership sets, we will prove that $ \tilde{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.

4.3.   Viability approach and a base numerical algorithm for computing leadership kernels

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 = (x_1, x_2, x_3, x_4) $ is the state vector and denote by $ f(x, u, v), \ v = (v_1, v_2, v_3), $ 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\to v(u) $ be a Borel measurable function with values in $ Q $. Consider the differential inclusion

Definition 4.1 (leadership property). A set $ K \subset R^n $ has leadership property if for any $ x_* \in K $ and any Borel measurable function $ v(\cdot) $ with values in $ Q $ there exists a solution $ x(\cdot) $ of (4.10) such that $ x(0) = x_* $ and $ x(t) \in K $ for all $ t \geq 0. $

Definition 4.2 (leadership kernel). For a given compact set $ G \subset R^n $ denote by $ Lead(G) $ the largest subset of $ G $ with the leadership property. This subset is called the leadership kernel of $ G $.

Let

be a family of state constraints. In the simulations below, $ g(x) = \max_ig_i(x) $ is set, where $ g_1(x) = |x_1-35|-7 $, $ g_2(x) = |x_2-47.5|-12.5 $, $ g_3(x) = |x_3|-0.2 $, $ g_4(x) = |x_4-5|-5, $ and $ g_5(x) = |\omega(x)| - \omega_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_{\lambda} $ is equivalent to the constraints (4.2) and (4.4) if $ \lambda = 0 $.

It is required to construct a function $ V $, such that

Denote by $ \delta > 0 $ the time step length. Let $ h: = (h_1, ..., h_n) $ be the space grid steps with $ |h| = \max\limits_{1\leq i \leq n}h_i $. For any continuous function $ \mathcal{V}: R^n \rightarrow R $ introduce the following upwind operator:

where $ f_i $ are the components of $ f $, and

Note that $ p_i^{right} $ and $ p_i^{left} $ are typical finite differences in Godunov-like methods for hyperbolic problems. The operator $ \Pi $ will be applied to grid functions to return grid functions.

Denote $ \mathcal{V}^h(x_{i_1}, ..., x_{i_n}) = \mathcal{V}(i_1h_1, ..., i_nh_n) $, $ g^h(x_{i_1}, ..., x_{i_n}) = g(i_1h_1, ..., i_nh_n) $ being the restrictions of $ \mathcal{V} $ and $ g $ to the grid. Let a sequence $ \{\delta_\ell\} $ is chosen, such that $ \delta_\ell \rightarrow 0 $, $ \sum_{\ell = 0}^{\infty}\delta_\ell = \infty $.

Consider the following grid scheme:

It can be proven that $ {\mathcal V}^h_\ell $ monotonically point-wise converges to a grid function $ V^{\delta, h} $, and this function approximates the function $ V $ defining leadership kernels according to formula (4.11) if $ \delta $ and $ h $ are sufficiently small.

Proposition 1 (see ^[30] for the proof). Let $ |{\mathcal F}_v(x)| \leq B $ for all $ x \in R^n $ and $ \delta_{\ell}/|h| \leq 1/(B\sqrt{n}) $ for all $ \ell. $ Then $ {\mathcal{V}}_{\ell}^h \nearrow {V}^h $ as $ \ell \rightarrow \infty $.

Note that in the case $ n = 4 $ the relation $ \delta/|h|\leq (2 B)^{-1} $ should hold. Other secondary stability conditions can be found in ^[31,21].

The estimate $ |{\mathcal{V}}^h - V| \leq C \sqrt{|h|} $ is expected (cf. ^[31,21]), and therefore $ {\mathcal{V}}_{\ell}^h $ approximates $ V $ if $ \ell $ is large and $ |h| $ is small. The stopping criterion in the computation is $ \sup\limits_{\substack{\hbox{ grid}}}|{\mathcal{V}}_{\ell+1}^h - {\mathcal{V}}_{\ell}^h| \leq \varepsilon $.

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 $ v_1, v_2 $, and $ v_3 $ 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:

where $ \bar{v}_i, i = 1, 2, 3$ , are new disturbances constrained in the same way as the old ones, cf. (4.3).

Second, put $ \tilde{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 $ x_1, x_2, x_3, x_4 $. The variables $ v_1, v_2 $, and $ v_3 $ are not considered as new states, but as accumulated optimal values of the new artificial controls.

4.4.   Modified (relaxed) algorithm for quick treatment of additional PT1 filters

Let $ w_i, i = 1, 2, 3$, be given grid functions. Consider now the following operator:

where

with the substitution

The recurrent computation scheme (4.12) transforms now into the following one:

Here, $ \bar{v}^{\ell, \;{ {max}}}_i(x), \;\; i = 1, 2, 3$ , are the maximizes, over $ \bar v_i, \;\; i = 1, 2, 3$ , when computing the term

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 $ w_i(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 $ \bar v_i, i = 1, 2, 3 $. Therefore, $ w^{\ell}_i(x), i = 1, 2, 3, $ accumulate, according to equation (4.13), values of $ \bar v_i, 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 $ w_i^0(x), i = 1, 2, 3, $ can be chosen arbitrary, under the constraints (4.3) with $ v_i $ being replaced by $ w^0_i(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 = \tilde{u}(x) $, for all admissible controls $ \bar{v}_i(t), i = 1, 2, 3 $. Therefore, level sets of the limiting function are leadership sets.

5.   Simulation results

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 $ v_1, v_2, v_3 $ will be checked.

First, the ability of the control $ \tilde{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 $ pCO_2 $ 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 = x_1(t) $, $ y = q\big(x_1(t), x_2(t), x_3(t), x_4(t)\big) $ and $ x = x_1(t) $, $ y = q\big(x_1(t), pCO_2^0, 0, p_v^*\big) $, 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 $ x_3 $ is shown. One can see a good closeness of $ x_3(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 $ x_3 $ (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\times 100\times 100 \times 100 $ cells was used, and about $ 10^4 $ 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 $ x_4 = p_v^* $ is shown. The vertical axis measures the discrepancy $ \omega $. 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.

6.   Discussion and Conclusion

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 $, $ pCO_2 $, 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 $ CO_2 $ 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 $ CO_2 $ 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]).

Acknowledgements

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.

Conflict of interest

The authors declare that there is no conflict of interest regarding the publication of this article.