Abstract
Impaired cerebral autoregulation leads to fluctuations in cerebral blood flow, which can be especially dangerous for immature brain of preterm newborns. In this paper, two mathematical models of cerebral autoregulation are discussed. The first one is an enhancement of a vascular model proposed by Piechnik et al. We extend this model by adding a polynomial dependence of the vascular radius on the arterial blood pressure and adjusting the polynomial coefficients to experimental data to gain the autoregulation behavior. Moreover, the inclusion of a Preisach hysteresis operator, simulating a hysteretic dependence of the cerebral blood flow on the arterial pressure, is tested. The second model couples the blood vessel system model by Piechnik et al. with an ordinary differential equation model of cerebral autoregulation by Ursino and Lodi. An optimal control setting is proposed for a simplified variant of this coupled model. The objective of the control is the maintenance of the autoregulatory function for a wider range of the arterial pressure. The control can be interpreted as the effect of a medicament changing the cerebral blood flow by, for example, dilation of blood vessels. Advanced numerical methods developed by the authors are applied for the numerical treatment of the control problem.
1. Introduction
Cerebral flow autoregulation is a process which aims to maintain proper and stable cerebral blood flow. By means of cerebral flow autoregulation, the body is able to deliver sufficient blood containing oxygen and nutrients to the brain tissue for metabolic needs and remove CO2 and other waste products. The most important objective of autoregulation is maintaining an appropriate level of brain perfusion, which is essential for life, since the brain has a high metabolic demand. It should be noticed that the brain is very sensitive to over- and underperfusion.
In recent years, mathematical modeling of autoregulation of cerebral blood flow (CBF) becomes very popular. Various methods like, for example, time- and frequency-domain analysis (see [1]) and principal dynamic modes (see [2]) are applied to analyze the dynamic relationship between the arterial blood pressure and cerebral blood flow velocity. Also, the so-called lumped parameter models based on the analogy to electric circuits (see, e.g., [3–5]) are widely used to simulate cerebral flow autoregulation mechanisms. Usually, the models are based on nonlinear ordinary differential equations (ODE) describing the flow balance in different parts of the cerebral vascular system (see, e.g., [6]). Autoregulation effects arise due to feedback relations governing the vascular volume (). It is supposed that is proportional to the square of mean vascular radius (), and the resistance () is inversely proportional to the fourth power of according to the Hagen-Poiseuille law. Therefore, , which yields the relation , or . It should be noticed that the experimental value of the exponent in the last relation is about (see [7]). This difference points out to the oversimplification of the cerebrovascular bed considered as the parallel arrangement of several equal microvessels. A more anatomically plausible, hierarchical model of the cerebrovascular bed is proposed in [7]. It consists of 19 compartments representing the whole range of vascular sizes and respective CO2 reactivities derived from literature data.
In this paper, we use the model from [7] in the context of autoregulation. We consider two methods for describing the dependence of the vascular radius on the arterial blood pressure. First, a polynomial dependence supported by fitting the corresponding coefficients to experimental data on CBF autoregulation in newborns is used. Introducing a hysteretic behavior into the resulting model is tested. Second, the ODE model from [6] is coupled with the model of [7]. A conflict control problem for such a model is considered. The control is interpreted as intake of a drug controlling the vascular volume, whereas unpredictable changes in arterial pressure are considered as disturbance.
Thus, the paper presented can be considered as a research demonstrating the idea to couple phenomenological compartment ODE models with physically consistent blood flow models derived from hydrodynamic laws. This opens the following perspectives. We are going to adapt the hydrodynamic part of the model to various types of cerebrovascular beds including germinal matrices of premature infants. This should be possible due to the flexibility of hydrodynamic models that can use precise descriptions of blood vessel networks, account for the assistance of vessel muscles in the blood circulation, or, conversely, assume muscular passivity of blood vessels. Moreover, different non-Newtonian and micropolar fluids can be used. The most important feature of hydrodynamic models is that they give the exact value of the blood pressure in all vessel types. Therefore, it becomes possible, based on the information about vessel sizes, to find the highest stress concentration in vessel walls. Then, using methods of elasticity and strength theory, breaking of the most fragile vessels can be predicted. On the other hand, control parameters included into the coupled model may be used to provide a regime that excludes damaging stresses in the most vulnerable vessels. Another peculiarity of the hydrodynamic part of the model is that only physical relationships between variables are imposed. Some artificial relations are not necessary. As for the phenomenological part of the model, related to ordinary differential equations, it should be refined, and phenomenological assumptions should be replaced by biophysical ones.
2. The Hierarchical Cerebrovascular Model
In this section, a hierarchical cerebrovascular model proposed in [7] is outlined. The derivation of expressions defining the CBF and the pressure in each compartment is explained. Notice that such a derivation is omitted in [7].
Let () be the number of arterial (venous) levels (see Figure 1), and . Denote by the pressure at the input of the th level. Notice that either (arteries) or (veins). Assume that and , where and are the arterial and venous pressures, respectively. Let , , , , and be the resistance, the number, the radius, the CO2 reactivity, and the length of vessels of the th level, respectively. The radius, , of vessels of the ith level depends on the partial CO2 pressure, , and is computed by the formula According to the Poiseuille law, the resistance of vessels is given by the formula , where is the dynamic viscosity of blood.
If CBF is known, the following relations obviously hold (Kirchhoff’s law): Therefore, CBF can be found from the following equality of pressures on the arteriovenous junction: Therefore, Finally, if we number all levels from the top down and remember that , we have Notice that, formally, and .
Combining formulae (1), (5), and (6) yields where are computed by formula (1).
Figure 2 shows the blood pressure depending on the vascular level. The computation is performed using formula (7). Here, we set , , and . The values of , , , and , , are taken from [7].
Figure 3 presents the flow velocity in each vessel depending on the level. The velocity is computed as , where CBF is given by formula (5). The condition is assumed.
Figure 4 shows the dependence of the pressure in the 8th, 9th, and 10th compartments on the partial CO2 pressure. The computation is performed using formula (7).
Figure 5 shows some statistical stability of the model with respect to the variation of the parameters , , , and . These parameters were randomly varied around the reference values reported in [7]. The uniform distribution of the variations was used, and the amplitude of the variations was equal to 20%. The dependence of the pressure on the 10th level on the partial CO2 pressure was computed in each test, and 200 random tests were done. The solid curve presents the average of the curves computed in each test. The dashed curve corresponds to the reference values of the parameters. The computations are done using formula (7).
3. Autoregulation Using a Polynomial Feedback
Now, the model described in the previous section is extended by adding a polynomial dependence of the vascular radiuses on the arterial blood pressure and adjusting the polynomial coefficients to experimental data to gain the autoregulation behavior. Additionally, a Preisach hysteresis operator is introduced into the resulting model to simulate the hysteretic dependence of the cerebral blood flow on the arterial pressure. Throughout this section, the assumption holds.
Let be a reference value of the arterial pressure and , , the vascular radiuses corresponding to . Assume that the dependence of on is chosen as follows: where is a parameter. Such a form of is motivated by the fact that they appear as fourth powers in the expression (5) so that the CBF becomes proportional to , and therefore the choice stabilizes CBF. Formally, this choice of yields the relation which implies (for all ) the equality To adapt the model to experimental data (see [8–10]), consider the following modification: Therefore, where is given by formula (11) and is a scale factor to adjust CBF to the age of infants. Set , , and , which corresponds to hemodynamic system of premature infants of 31–34 weeks’ gestational age with 260 g brain weight and of 15.5 mL/100 g/min (cf. [9, 10]). The values of , , and , , , are taken from [7]. The coefficients , and are fitted through the minimization of the residual where the pairs , are chosen according to [9] as follows: The values (in [1/Pa]) of the minimizers of (13) read It should be noticed that the polynomial appearing in square brackets in (11) with coefficients (15) represents a sigmoidal function providing autoregulation. Thus, the polynomial ansatz is a method of constructing an appropriate sigmoidal function.
Figure 6 shows the effect of autoregulation when using the feedback law given by formulae (11) and (12) with the coefficients (15).
Figure 7 shows the same simulation as in Figure 6, but a hysteretic behavior of the process is added by applying a Preisach hysteresis operator. Thus, if the arterial pressure increases, the output follows the main path (cf. Figure 6). If the pressure starts to go down at some point, the output path goes back below the main path. If the pressure again increases, the output returns to the main path and follows it. Notice that Preisach hysteretic relations express wetting/drainage processes in structures containing a large number of small vessels and capillaries. A fast numerical method implementing the Preisach hysteresis operator is developed in [11]. The simulation of hysteretic behavior was motivated by experimental results reported in papers [12–14].
In [12], an asymmetry of autoregulation during spontaneous increases and decreases of cerebral perfusion pressure (CPP) in brain injury patients was experimentally shown. The autoregulatory response was significantly greater during increase than during decrease in CPP.
Experimental data of [13] based on a pharmacological approach to evaluate dynamic cerebral autoregulation (dCA) gain to transient hypotension and hypertension in healthy patients showed that cerebral autoregulation is different for rising and falling blood pressure. Namely, dCA gain to transient hypotension was consistently greater than dCA gain to transient hypertension.
According to [14], a strongly asymmetric dynamic response of the cerebral autoregulation was seen in the majority of patients with head injury. It might also have been present, albeit to a lesser degree, in the normal subjects. The findings suggest that nonlinear effects may be present in the operation of the cerebral autoregulation mechanism.
4. Autoregulation Using an ODE Model by Ursino and Lodi
In this section, the ODE model from [6] is coupled with the model of [7] (cf. Section 2). The aim of this section is to compare simulation results produced by our coupled model with the results of [6]. Thus, the data reported in [6] are used in our coupled model. These data correspond to adult persons, and therefore the autoregulation area approximately extends from 50 to 140 [mmHg], and CBF is around 11.4 [mL/s].
The variables and parameters used in [6] are the following: : the arterial pressure (input), : the dural sinus pressure (parameter) = 6.0 mmHg, : the partial pressure of CO2 (parameter) = 0, 10, 20, 30, and 40 mmHg, : arterial-arteriolar volume (variable), : the intracranial pressure (variable), : capillary pressure (variable), : the arterial compliance (variable), : the arterial resistance (variable), : the cerebral blood flow (CBF) (variable), : the CBF deviation normalized to the basal level (variable), : a sigmoidal static function, : the amplitude of the sigmoidal curve, see (24), : the central slope = , see (24), : the maximum autoregulation gain = 1.5 mL/mmHg, : 0.75 mL/mmHg, : 0.075 mL/mmHg, : the central value of the sigmoidal curve = 0.15 mL/mmHg, : the resistance to the formation of cerebrospinal fluid = mmHg·s/mL, : the proximal venous resistance mmHg·s/mL, : the resistance to the outflow of cerebrospinal fluid = 526.3 mmHg·s/mL, : the injection rate of cerebrospinal fluid = 1/30 mL/s, : the value of CBF required by tissue metabolism = 12.5 mL/s, : the time constant of the regulation = 20 s, : a fitting constant = 0.2, : 0.11 mL−1.The new model consists of the same ODEs as the model from [6], and the relations Here, the relations (19) and (22) replace oversimplified formulas for and from [6]. The definitions of the functions CBF (returns the cerebral blood flow) and (returns the vascular resistance) are taken from Section 2. They are given by the formulas where , , , and are the number, the length, the radius, and the CO2 reactivities of vessels of the th level, respectively (see [7]).
The function is defined by the relations The initial state is mL/mmHg and mmHg.
Thus, the difference with the model developed in [6] consists in the computation of the cerebral blood flow and the arterial resistance using the model of cerebrovascular bed considered in [7]. In contrast to [6], and are defined by the relations (19) and (22), respectively. Additionally, it was assumed that the vascular radiuses are proportional to the square root of the arterial-arteriolar volume; that is, in (23).
Figure 8 shows two outputs of the model (16) and (17) versus . The input signal is a linear function with the slope equal to 2/3 mmHg/s, and the partial pressure of CO2 equals zero.
Figure 9 presents five outputs of the model (16) and (17) corresponding to the following magnitudes of the partial CO2 pressure: 0, 10, 20, 30, and 40 [mmHg]. The input signal is the same as in the previous simulations.
5. Conflict Control Setting
The crucial variable defining the behavior of the system (16)–(22) is the arterial-arteriolar volume controlled by the arterial compliance ; see (20). The behavior of can be influenced by a control parameter, , added to the right-hand-side of (16). The effect of can be interpreted as the intake of a drug increasing/decreasing the response of the CBF autoregulation system to the deviation of . Moreover, assume now that the realization of is formed by an interfering factor that intends to crash the nominal . Assume that this factor “chooses” the rate, , of at each time instant. Additionally, observe that simulations of the model (16)–(22) show small variations of the variable , and therefore (17) can be preliminary neglected. Summarizing these assumptions, arrive at the following conflict control system: where the function is defined by formulae (18), (19), and (20) with mmHg. We assume that the control, , and the disturbance, , are restricted as follows: and . The objective of is the minimization of the functional whereas strives to maximize it. Here, is the termination time of the control process, and is a cutoff function such that if and it rapidly falls to zero outside of the above interval. Notice that the usage of means that the deviation is not assumed to be small outside of the reasonable range of the arterial pressure.
The system (25) contains nonlinear dependencies given by a computer subroutine (cf. the computation of ). Moreover, a disturbance is involved into the dynamics. Thus, the application of traditional control design methods based on Pontryagin’s maximum principle is not applicable. Nevertheless, the dynamical programming principle related to Hamilton-Jacobi-Bellman-Isaacs (HJBI) equations is suitable in this case. The application of this technique requires stable grid methods for solving HJBI equations arising from conflict control problems. Such methods are developed by the authors (see [15, 16]) and can be applied to the problem (25)-(26). The next section outlines an upwind grid method for the solution of HJBI equations.
5.1. Finite-Difference Scheme
Notice that the state vector of the problem (25)-(26) is two-dimensional. For generality, consider -dimensional conflict control problem with the dynamics and the gain defined as follows: Here, and are the control parameters of the minimizing and maximizing players, respectively, and are feedback strategies of the players, and is the start position of the control process.
Let and be time and space discretization step sizes. Let be a grid function. Define an upwind operator, , as follows: where and runs over all grid points.
Let . Denote , , and introduce the following notation: that is, runs over all grid points. The numerical scheme is the following:
The proof of the convergence of to the gain is given in [15, 16]. Optimal feedback strategies of the players can be computed by storing a minimizer and a maximizer in (29) at each sampled time instant and each grid point. Another way of the construction of optimal feedback strategies is the so-called extremal aiming procedure (see, e.g., [16] for an explanation).
5.2. Result
Figure 10 shows the result of the application of the above sketched technique to the minimization of the functional (26). The input signal grew from zero to 150 [mmHg] with a constant rate and then fell to 80 [mmHg]. One can observe the extension of the horizontal plateau in the case of applying the control. Thus, the computed feedback control, , where is measured and is defined from the model (25), allows us to fit the intake of the drug to pressure jumps.
6. Conclusions
Simple models based on Kirchhoff’s law and Hagen-Poiseuille flow of Newtonian fluids are studied here in the context of CBF autoregulation. Such models do not require a great amount of data. Information on the approximate number of vessels, their size, length, and reactivity is already available for cerebral vascular system including germinal matrices of preterm newborns. Moreover, the model behavior is stable with respect to the variation of these data. Using such models allows us to estimate pressures in the germinal matrix and their dependency on the partial CO2 pressure. Coupling this model with models of CBF autoregulation may be useful to gain a better understanding of mechanisms of violation of hemodynamics in germinal matrices.
The future work will be related to the enhancement of the above described coupled model. First, a method for generation of hierarchical cerebrovascular networks will be developed. Then, the hydrodynamic part of the model will be extended by accounting for non-Newtonian and micropolar properties of blood (modification of Poiseuille’s law). Expansion/contraction travelling waves propagating along vessels and the curvature of vessels will be included in the model. Such an enhanced model should predict the precise pressure distribution over the vessel network. The information about the length, radius, and wall thickness of vessels will allow us to find vessels having maximal stress values in their walls. These results will be used in finite element simulations to estimate the damage risk of such vessels. The optimal control setting will be enhanced in such a way that the objective will consist in keeping stresses in a safe range over all vessels.
Conflict of Interests
The authors declare that they have no conflict of interests.
Authors’ Contribution
Renée Lampe, Nikolai Botkin, and Varvara Turova equally contributed to this work.
Acknowledgments
The authors would like to thank two anonymous referees for their valuable comments and suggestions. Financial support of the Buhl-Strohmaier-Stiftung and Stiftung Würth is gratefully acknowledged.