Optimal control problem arising in mathematical modeling of cerebral vascular pathology embolization

Arteriovenous malformation (AVM) of the brain is a congenital vascular abnormality, in which the arterial and venous blood pools are intertwined and directly connected. This dangerous disease causes a high risk of intracranial hemorrhage and disrupts brain functioning. The preferred method of AVM treating is embolization, which is the endovascular filling of abnormal AVM vessels with a special embolic agent. Despite the fact that this method is widely used in neurosurgery, in some cases its use is accompanied by perioperative AVM vessels rupture. In this regard, the aim of this work is to study the optimal scenarios for multi-stage AVM embolization from the effectiveness and safety of the procedure point of view. Mathematically, the joint movement of blood and embolic agent in the AVM body is described on the basis of a one-dimensional two-phase filtration model, which takes into account the redistribution of blood to surrounding healthy vessels. For the numerical solution of the resulting integro-differential system of equations, a monotonic modification of the CABARET scheme is used. To find optimal embolization scenarios, the optimal control problem with phase constraints arising from medicine is formulated. A modified particle swarm optimization method is used to solve this problem numerically. This technique is used to obtain optimal embolization scenarios on the basis of real patients clinical data collected during neurosurgical operations.


Optimal control problem arising in mathematical modeling of cerebral vascular pathology embolization Tatiana Sharifullina * , Alexander Cherevko & Vladimir Ostapenko
Arteriovenous malformation (AVM) of the brain is a congenital vascular abnormality, in which the arterial and venous blood pools are intertwined and directly connected. This dangerous disease causes a high risk of intracranial hemorrhage and disrupts brain functioning. The preferred method of AVM treating is embolization, which is the endovascular filling of abnormal AVM vessels with a special embolic agent. Despite the fact that this method is widely used in neurosurgery, in some cases its use is accompanied by perioperative AVM vessels rupture. In this regard, the aim of this work is to study the optimal scenarios for multi-stage AVM embolization from the effectiveness and safety of the procedure point of view. Mathematically, the joint movement of blood and embolic agent in the AVM body is described on the basis of a one-dimensional two-phase filtration model, which takes into account the redistribution of blood to surrounding healthy vessels. For the numerical solution of the resulting integro-differential system of equations, a monotonic modification of the CABARET scheme is used. To find optimal embolization scenarios, the optimal control problem with phase constraints arising from medicine is formulated. A modified particle swarm optimization method is used to solve this problem numerically. This technique is used to obtain optimal embolization scenarios on the basis of real patients clinical data collected during neurosurgical operations.
Arteriovenous malformation (AVM) of the brain is a congenital vascular abnormality, in which the arterial and venous blood pools are intertwined and directly connected, bypassing the capillary vessel network. By type of pathological blood vessels, connecting the arterial and venous pools, AVMs can be divided into racemic pathology (consisting of large number of small diameter vessels, which are chaotically intertwining and intersecting with each other) and fistula one (consisting of large vessels). The resistance of the corresponding part of circulatory system is reduced due to the capillary network absence. This causes changes in both the hemodynamic parameters (flow rate and pressure) and the strength properties of the blood vessels [1][2][3][4] .
Currently, the most common methods of AVM treatment are microsurgical, endovascular and radiosurgical 5 . AVMs have a complex effect on intracranial hemodynamics, which varies with changes in morphology and angioarchitecture, especially with different treatment methods. Endovascular embolization is the filling of abnormal blood vessels with a special embolic agent in order to shut them off from the bloodstream. Due to the minimal invasiveness and the possibility of operating in deep, functionally significant brain areas, embolization is one of the most effective methods at the current level of medical development 6 . Despite the fact that this method is widely used, the risk of perioperative vascular rupture is still a serious danger 7 . Therefore, the study of optimal scenarios for multi-stage AVM embolization from the effectiveness and safety of neurosurgical surgery point of view remains an urgent task.
For the mathematical description of blood flow, a one-dimensional approximation of the Navier-Stokes equations is widely used, obtained by averaging these equations over blood vessel cross-section 8 . There are works in which the equations describing flows in pipes with rigid and elastic boundaries are analytically investigated [8][9][10] . Viscous fluid flow in a network of soft tubes is also modeled on the basis of a one-dimensional approximation of mass and momentum conservation laws [11][12][13] . Numerical simulation of hemodynamics for large blood vessels based on 3D-1D coupled flow is also considered [14][15][16] .
The interaction of AVM and blood flow in the surrounding vessels is often studied in analogy with electric and hydraulic networks [17][18][19] . Such models allow us to assess the impact of various embolization scenarios on Scientific Reports | (2022) 12:1302 | https://doi.org/10.1038/s41598-022-05231-w www.nature.com/scientificreports/ blood flow redistribution and are consistent with the general medical view on pathology hemodynamics. Various approaches are used to describe the joint flow of blood and embolic agent, for example, a two-phase flow model was used to simulate a viscous fluid drop movement of embolic agent through a bifurcation point or multifurcation point based on dimensionless Navier-Stokes equations for incompressible liquids 20 . A model of AVM embolization was proposed based on the concept of two-fluid modeling and scalar transport, where embolic agent and blood interaction along with its hardening is simulated by viscosity increase 21 . Another approach to the description of the embolization process is based on the two-phase filtration model in the authors' work 22 . The present paper is a extension of this work. In comparison to the present work in 22 blood flow into the adjacent healthy vessels was not considered, the cross-section area and porosity were assumed to be constant, and embolization was modeled without taking into account embolic agent solidification between the stages of surgery (simplified one-stage case). The optimal control problem 22 differs in form from the one described below, but is essentially the same, namely, it requires performing the operation in the most efficient way in compliance with safety constraints. The formulation of the optimal control problem in present paper is obtained from the control problem 22 by interchanging the constraint and the objective functional, which is typical of extremum search problems with constraints. In the work 22 , the search for an extremum of the functional was based on the calculation of its value on the grid in the three-dimensional control parameter space with subsequent interpolation. While present work uses the global optimization method -the particle swarm method, which significantly accelerates calculations, increases their accuracy and allows us to explore the control parameter space of a larger dimension. Despite the large number of publications dedicated to the AVM hemodynamics, research expansion in this direction is a great scientific interest. In this paper, mathematical modeling of AVM embolization is carried out as follows: the flow of blood and embolic agent through the AVM body is considered as filtration flow based on Darcy's law, taking into account the redistribution of blood to healthy vessels. This approach is justified for small-vascular racemic AVM compartments embolization description. It is assumed that both liquids are Newtonian with constant viscosity, incompressible and immiscible, while ignoring the interfacial capillary forces, the effect of embolic agent adsorption, as well as the deformation and permeability of the vessel walls.
As a result, the mathematical description of the joint movement of blood and embolic agent in the AVM body is based on a one-dimensional two-phase filtration model that takes into account the redistribution of blood to surrounding healthy vessels. For the numerical solution of the resulting system of integro-differential equations, a monotonic modification of the CABARET scheme is used, which localizes with high accuracy the strong and weak discontinuities that occur when solving this problem 23,24 . For numerical calculations real blood and embolic agent (ONYX18) viscosities were used. To find optimal embolization scenarios, the optimal control problem of multi-stage embolization with phase constraints arising from medicine is formulated. To solve this problem numerically a particle swarm optimization method is used 25,26 , which is specially modified in order to adapt it to the problem under consideration. The proposed technique is used to obtain optimal embolization scenarios on the basis of real patients clinical data collected during neurosurgical operations 22,27 .
In the next section, a mathematical formulation of the one-stage embolization problem is given, leading to an initial-boundary value problem for a system of integro-differential equations, including a hyperbolic partial differential equation, which is a quasilinear scalar conservation law with a nonconvex flux. Then the optimal control problem of one-stage AVM embolization is formulated, where the control is a time function setting the boundary condition of the problem and is included in the coefficients of the equation. The objective functional is the integral of the initial-boundary value problem solution at terminal time. The phase constraints on the control function represent additional integral and algebraic relations. The objective functional and constraints are selected according to medical indications. We introduce a class of piecewise linear continuous and discontinuous functions used as controls for the AVM embolization process. Further, the optimal control problem of multi-stage AVM embolization is formulated, which leads to a chain of integro-differential initial-boundary value problems, in which each subsequent problem uses the results of solving the previous problem as parameters of a new system of equations.
Next, we present a monotonic modification of the CABARET scheme, which localizes with high accuracy the strong and weak discontinuities that arise when solving the initial-boundary value problem under consideration. The modification of the particle swarm optimization method is proposed that is specially adapted to the considering problem in order to more fully investigate the admissible parameters area and accelerate the progress of the particle swarm to the absolute minimum. Numerical verification of the optimal control problem of multi-stage embolization using real patients clinical data was performed, on the basis of which AVM embolization optimal scenarios were obtained. Also we discusses the limitations of the proposed method, as well as the possibilities for its further development.

One-stage embolization problem formulation
In this paper, blood and embolic agent movement through small-vascular racemic AVM is modeled (Fig. 1) as a two-phase filtration process of immiscible incompressible liquids, where the displaced phase is blood, and the displacing phase is embolic agent. We also take into account the redistribution of blood to surrounding healthy vessels. We will consider a one-dimensional case, which mathematical description is based on Darcy's law and the mass conservation laws for each phase. In the one-stage embolization case, an initial-boundary value problem is obtained for a system of integro-differential equations, in which the hyperbolic partial differential equation is a quasi-linear scalar conservation law with a nonconvex flux function. This conservation law is similar to the one widely used to describe the multiphase filtration process 28,29 .
Based on this initial-boundary value problem, it is set the task of finding the optimal regime of one-stage AVM embolization from the effectiveness and safety of neurosurgical surgery point of view. The objective functional www.nature.com/scientificreports/ and constraints that arise in the optimal control problem are selected according to medical indications. The control is a time function that determines the embolic agent injection at the AVM input. The embolization effectiveness is determined by the degree of blood displacement from the AVM by the embolic agent at the operation end. One of the possible risks of the operation complications is the embolic agent entering in the vein, as a which result the AVM or vein vessels rupture may occur. Also, during the operation, vascular network has an additional stress, since when the embolic agent is injected, the pressure in the arterial vessels increases. Excessive pressure build-up is dangerous. Both of these factors increase the brain hemorrhage risk 30 and will be further used in the optimal embolization control problem.
Derivation of the conservation law simulating the process of one-stage embolization. Consider the one-dimensional joined filtration flow of blood and embolic agent through the AVM body. Denote by A the AVM cross-section area available for the phase flows, by m the AVM body porosity, by Q the volume flow rate of two-phase mixture. These values can vary both in space and in time. Let S(t, x) ∈ [0, 1] be the local saturation (blood concentration) of the porous medium by the displaced phase (i.e. the fraction of the pore volume occupied by blood). Taking this into account, the concentration of the displacing phase is 1 − S(t, x) . We introduce the Buckley-Leverett saturation function f(S) 29 , which sets the distribution of the phases volume flow inside a porous medium and is equal to the flow fraction of displaced phase in the total two phases flow. We derive the differential form of the integral conservation law of blood volume, assuming that the AVM walls are impermeable and the phases are incompressible fluids. The change in blood volume V 1 on the segment This change is due to the blood flow through the cross-sections at points x = x 1 and x = x 2 : Equating the values (1) and (2), we obtain the integral conservation law of blood volume where m(x) and A(x) are given piecewise continuous functions of the spatial coordinate x. The differential form of the conservation law (3) has the form Since the Eq. (4) is hyperbolic, its solution S(t, x) is considered in the class of piecewise continuous functions. We will assume that the Buckley-Leverett function f (S) ∈ C 2 [0, 1] satisfies the following conditions where ξ is the inflection point. Since the AVM side surface is impermeable and the liquids are incompressible, the total volume flow rate of the two phases Q(t) does not depend on the spatial coordinate x and is determined by the relation where Q b = Q f (S) is the volume blood flow and Q e = Q (1 − f (S)) is the volume flow rate of an embolic agent inside AVM. With this assumptions, the Eq. (4) can be written in the following equivalent non-divergent form (3) ; and Q e -volume flow rates of blood and embolic agent inside AVM; q b and q e -volume flow rates of blood and embolic agent at the AVM input; p 1 -blood pressure at the AVM input; p 2 -blood pressure at the AVM output. Initial and boundary conditions. The differential equation (5) is solved in the domain where L and T are some positive numbers. Since at the initial time t = 0 there is only blood inside the AVM, then the initial value of its concentration is given by the following formulas Blood flow changes during the cardiac cycle is not taken into account due to the duration of the embolization process is significantly longer than the cardiac cycle duration, and the blood flow rate q b (t) = Q b (t, 0) at the AVM input is determined by the formula in which q b (S) is a given function that determines the redistribution of blood between AVM and surrounding healthy vessels, where is the average blood concentration in the AVM. At the AVM input ( x = 0 ) the embolic agent flow q e (t) = Q e (t, 0) determines the embolization process and is a control function. Then the flow rate Q(t) of two-phase mixture inside the AVM, included in (5), is determined by the formula Assuming that the blood and embolic agent flows satisfy the inequalities we obtain that the flow of two-phase mixture in consequence, the characteristics of the differential equation (5) propagate with velocities Q(t)f ′ (S)/(m(x) A(x)) > 0 . This means that for the correct formulation of the initial-boundary value problem for the Eq. (5) one boundary condition must be set on the left boundary of the segment [0, L] and there is no need to set boundary conditions on the right boundary of this segment. It follows from formula f (S(t, x)) = Q b (t, x)/Q(t) and relation (8) that the boundary value S(t, 0) is expressed in terms of boundary phase flows and is defined as follows where the function S (t) is given by the integral formula (7). As a result, we obtain the integro-differential problem (5), (6), (8), (9) in which the integral relation (9) plays the role of a non-standard boundary condition for the function S(t, x). It is important to emphasize that the control function q e (t) is included both in the coefficients of the Eq. (5) and in the boundary condition (9).

Optimal control problem in the case of one-stage embolization
To formulate the embolization optimal control problem, it is necessary to determine objective functional and constraints arising from medical indications, which establish the effectiveness criteria and safety conditions of the operation. The embolization effectiveness will be determined by how completely the embolic agent fills the AVM at the operation end. Mathematically, this condition lead to minimizing the objective functional where S q e (t, x) is the solution of the initial-boundary value problem (5), (6), (8), (9) with control function q e (t) , which sets the embolic agent flow at the AVM input.
The safety embolization conditions can be chosen in the form of the following restrictions , www.nature.com/scientificreports/ where p 1 is the pressure at the AVM input, and p * is the specified critical pressure. Fulfillment the condition (11) avoids the penetration of the embolic agent into the venous bed, and fulfillment the inequality (12) avoids an excessive increase in blood pressure in the AVM. Violation of the conditions (11) and (12) increases the cerebral hemorrhage risk 30 .
Taking into account the fact that Q b = Q f (S) , the pressure function is determined from Darcy's law for AVM 29 where p 2 (S) is a given function of the blood pressure at the AVM output. The local resistance to blood flow inside the AVM is given by the ratio where η b is the blood viscosity, K is the absolute permeability of the AVM and k b (S) is the blood relative phase permeability; the values of η b , K and k b (S) are assumed to be known. For the given functions f(S), r b (S) , q b (S) and p 2 (S) , Eqs. (5), (8) and (13) form a closed system of integro-differential equations for determining the blood concentration S(t, x) inside the AVM, the mixture flow rate Q(t) and the pressure p 1 (t) at the AVM input. The pressure p 1 (t) is used to verify the considered optimal AVM embolization model by comparing with the pressure obtained during intraoperative monitoring. Thus, in the optimal control problem, it is necessary to find a function q e (t) that minimizes the objective functional (10), which is determined on the basis of the problem (5), (6), (8), (9) solution. Besides, it must be done the restrictions (11), (12) on the blood concentration S and the pressure p 1 (t) determined by the formula (13).

Control functions class
We consider a numerical solution of the problem (5), (8), (13) with initial (6) and boundary (9) conditions that satisfies the constraints (11) and (12). The choice of the control class is determined by one of the methods of embolization in real surgical practice: at the initial time interval, there is increase in the embolic agent supply; at the next time interval, a constant value of this flow is maintained; and at the last short time interval, there is a decrease in the embolic agent supply to the AVM input. We use linear interpolation of this embolization method. As a result, we consider the following continuous functions as a control function class where γ > 0 is a dimensionless parameter, Q(0) = q b (0) is the blood flow through the AVM before surgery, is a dimensionless function in which the parameter θ satisfies the condition 0 < θ ≤ T − ε , and ε ∈ (0, T) is a small constant. Admissible control class is expanded to include the function that is discontinuous at the initial moment of time obtained from the function (15) when θ → 0.
For the control class (14) and (15), time interval [0, θ] corresponds to an increase in the supply of embolic agent to the value γ Q(0) ; time interval [θ, T − ε] corresponds to maintaining a constant value of this flow; time interval [T − ε, T] corresponds to a decrease in the embolic agent supply to the AVM input. The control class (14) and (16) corresponds to the case when there is a discontinuous increase in the embolic agent supply.
As a result, the optimal embolization problem is led to finding the control parameters that determine embolic agent supply (14). The corresponding solution of the initial-boundary value problem (5), (6), (8), (9), (13) gives an absolute minimum to the functional (10), which corresponds to minimizing the average blood concentration in the AVM at the operation end. The solutions to the optimal control problem must satisfy the constraints (11) and (12), where (11) prohibits the embolic agent from entering the vein, and (12) restricts the blood pressure at the entrance to the AVM during embolization. Next, instead of the parameters (17), we will use the following control parameters www.nature.com/scientificreports/ where the parameter θ 1 corresponds to the time of increasing the embolic agent supply to the AVM input, θ 2 corresponds to the time of maintaining the achieved flow of embolic agent, γ is the intensity of the embolic agent supply.

Multi-stage embolization
Multi-stage embolization problem formulation. A real neurosurgical operation for AVM embolization usually proceeds in several stages, the number of which is denoted by N. Mathematical modeling of multistage embolization is a generalization of the one-stage embolization described in the previous section. Between the stages there are time breaks, during which the embolic agent remains fixes inside the AVM and does not change its location in the following stages of embolization. For the mathematical modeling, we assume that the duration of these time breaks is zero. Taking this into account, embolization at each i-th stage occurs on the time interval [T i−1 , T i ] , where T 0 = 0 . The cross-section area available for the flow of phases at the i-th stage is denoted by A i (x) , so A 1 = A is the original cross-section area. The fraction of cross-section A i occupied by blood is denoted by S i (t, x) , and the fraction of the original cross-section A 1 occupied by blood is denoted by Between the stages there is a change in the fraction of the AVM cross section available for the flow of phases which establishes a connection between the successive embolization stages. We will assume that in the AVM part available for the phases flow it is not change from stage to stage the filtration characteristics: porosity m, absolute permeability K, relative phase permeability k b and Buckley-Leverett function f. In the multi-stage case, the embolization process at each i-th stage is modeled similarly to the one-stage embolization case, which was described above. Thus, the initial-boundary value problem simulating the i-th embolization stage on the time interval [T i−1 , T i ] , is similar to (5) For multi-stage embolization, we consider a control that at each stage similar to the control for one-stage embolization (14)- (16). Such control can be considered as a linear approximation of the clinical process of multi-stage embolization and written as follows . obtained from functions (28) when θ i → 0. As a result, the optimal embolization problem is reduced to finding the control parameters that determine the embolic agent supply (27), and the corresponding solution of the initial-boundary value problem of multi-stage embolization (18)- (23) gives an absolute minimum to the functional (24), which corresponds to minimizing the average blood concentration in the AVM at the multi-stage operation end. The solutions to the multi-stage optimal control problem must satisfy the constraints (25) and (26), where (25) prohibits the embolic agent from entering the vein, and (26) restricts the blood pressure at the entrance to the AVM at i-th embolization stage. Next, instead of the parameters (30) we will use the following control parameters where the parameter θ 1i corresponds to the time of increasing the embolic agent supply to the AVM input, θ 2i corresponds to the time of maintaining the achieved flow of embolic agent, γ i is the intensity of the embolic agent supply at the i-th embolization stage.

Methods
Finite difference scheme. For the numerical solution of the initial boundary value problem (22), (23) for the system of integro-differential equations (18)-(21), we will use a monotone modification of the explicit two-layer time the CABARET scheme 31 , which has a second-order approximation on smooth solutions. This modification of the CABARET scheme was proposed in paper 23 and is used to model the problems of hemodynamics 22,32 . The main advantages of this scheme are due to the fact that it is given on a compact spatial stencil and, when approximating linear equations, is time-reversible and accurate with two different Courant numbers r = 0.5, 1 , which gives it unique dissipative and dispersive properties 33 .
Numerical simulation is carried out one after another at all stages of embolization, and the same difference scheme is applied at each stages. We present the algorithm of CABARET scheme for the i-th stage. Let us divide the segment [0, L] on which the considered problem is solved into J equal parts of length h = L/J and define a rectangular uniform difference grid where h is a constant spatial step, and τ n i is the time step at the n-th time layer, determined from the CFLstability condition at the i-th embolization stage. Let us denote by S c i the approximate numerical value of the blood concentration S i at the i-th stage. The CABARET scheme uses the flux (u i ) n j = S c i (t n i , x j ) and conservative (U i ) n j+1/2 = S c i (t n i , x j+1/2 ) variables set respectively in the integer x j and half-integer x j+1/2 = x j + h/2 spatial nodes of the difference grid (32).
Let   The multi-stage embolization optimal regime: modified particle swarm optimization method. To solve the problem of optimal multi-stage embolization, we use a numerical algorithm based on the particle swarm method 25 . To apply this method, the optimal control problem (24)-(26) is rewritten as a functional with penalty functionals I 1 and I 2 , where I 1 [S i ] = 0 when the constraint (25) is met and I 1 [S i ] = 1 when the constraint is violated; similarly, I 2 [S i ] = 0 when constraint (26) is met and I 2 [S i ] = 1 when this constraint is violated; here the solutions S i depends on the control function q e (t) . In the following calculations, in formula (40) we assume R = 1000 . In the future, the objective functional modification (40) makes it possible to significantly simplify the application of the particle swarm method near a priori unknown part of the boundary domain of control parameters set by the constraints (25) and (26). The required parameters θ 1i , θ 2i and γ i , where i ∈ 1, N , included in the functional (40), are represented as the coordinates of the abstract particle x ∈ R 3N , which set M = {x 1 , . . . , x M } is called as the particles swarm 25,26 . The particle swarm method for minimizing the functional consists in setting the initial position of each particle and organizing the iterative process of moving these particles so that after a given finite number of iterations I all the particles end up in a small neighborhood of some point x * that is the global minimum to the functional (40). At each iteration it is determined the displacement velocity v j and the best position p j of each particle x j , and also the global best position p g of all the particles for all the previous iterations. At the end of the iterative process, point x * = p g is the solution of the minimization problem.
The search for the vector x * is as follows. First, in the domain of control parameters (31) we allocate M bounded disjoint subdomains j ; for each j the initial position of single particle is set according to the formula is a functional that randomly selects the point in the domain j . The initial velocity of each particle is v 0 j = 0 and its initial best position is p 0 j = x 0 j . For all particles x 0 j the value of the functional Ĵ N (x 0 j ) is calculated, after which the formula www.nature.com/scientificreports/ is used to find the initial global best position for all particles.
Let the values x k j , v k j , p k j , j = 1, M , and p k g be known at the k-th step of the iteration. The values at the (k + 1) -th step of the iteration are determined in several stages. When k + 1 = I the iterative process is ended. The proposed algorithm differs from the classical particle swarm optimization method 25 in that a special initial particle distribution is selected and that at each algorithm iteration the positions of all particles are reinitialized with a certain probability P. This modification was required for a more complete examination of the admissible parameter domain and to accelerate the progress of the swarm to an absolute minimum.
Using clinical data of real patients. To verify the optimal control embolization problem, we use clinical data obtained during monitoring of hemodynamic parameters in process of neurosurgical operations. In these operations, total embolization of the AVM was achieved, that is, the AVM available volume at the operation end was completely filled with the embolic agent. The monitoring was carried out at the National Medical Research Center named after academic Meshalkin 27 using a Philips ComboMap system and the Philips ComboWire sensor (sensor diameter is 0.36 mm and its length is 1.85 m), which measures blood velocity and pressure inside the cerebral vessels near the pathology. This way the pressure and velocity at the AVM input were obtained before, during and after the embolization, as well as at the AVM output before and after the operation. To obtain the geometrical AVM parameters: the length L, the cross-section area A as well as input artery cross-section area, data from angiographic studies (perioperative X-ray tomography) were used (Table 1). In numerical calculations, the blood viscosity is η b = 4 cP , and the embolic agent viscosity is η e = 18 cP , which corresponds to the viscosity of the common embolic agent ONYX18 6 , where cP is viscosity unit centipoise.
Using the method proposed in the previous papers 22 , the functions k b (S) , p 2 (S) , q b (S) , f(S) and r b (S) are constructed for each patient on the basis of clinical data. These functions are included in the initial-boundary value problem for a system of integro-differential equations (18)- (23). In the numerical solution of the optimal control problem (24)-(26), the value p * is chosen as the maximum pressure at the AVM input during the operation.
In To verify the proposed model, we consider the multi-stage embolization problems for patients K, S and P. The geometrical and filtration characteristics of the pathologies of these patients used for calculations are shown in Table 1 and Fig. 2.
The investigation was conducted in accordance with the Helsinki Declaration and approved by the Inspection Commission of the Meshalkin Clinic (Session No. 549). All included patients gave their informed consent.

Results
Numerical calculations were performed for patients K, S and P based on the clinical data given above. For the calculations, a space-uniform difference grid was used (32), where J = 100 , with a time step given by the formula (33), in which CFL number r = 0.5 . The small parameter included in the control (28), (29) is ε = 0.1 s, which corresponds to the medical practice of quickly stopping the embolic agent supply at the end of each embolization stage. According to the papers 25,26 , the parameters of the particle swarm algorithm were chosen as follows: P = 0.15, w 1 = 0.9, w 2 = 0.5, c 1 = c 2 = 1.49 . At each stage, the number of iterations of the particle swarm method is I = 300.
In order to be consistent with the clinical data, numerical calculations were performed in such a way that the objective functional value was Ĵ N (T N ) ≤ 0.05 . It was shown that such a functional value for patients K, S is achieved in two embolization stages, and for patient P in three stages. A further increase in the number of embolization stages leads to smaller values of the objective functional. It should be noted that in a real neurosurgical operation, the number of stages can be more than three. To solve the optimal control problem of two-stage embolization by the particle swarm optimization method, the number of particles is chosen to be M = 16 , and for the problem of three-stage embolization, the number of particles is increased to M = 32 . For the optimal regimes obtained as a result of the numerical solution of the multi-stage embolization problem, Table 2 shows the values of the objective functional Ĵ N and the corresponding values of the parameters p g for three patients K, S and P. Figure 3 shows a comparison of the clinical and calculated pressures at the AVM input for the optimal regimes. This figure shows a good correlation between the calculated and clinical pressures for patients K and S. At the same time, for the patient P, the behavior of the clinical and calculated pressures differs to a greater extent. This probably due to the fact that this patient has a difficult clinical case, as indicated by the rupture of AVM vessels after surgery. Thus, it can be assumed that the process of clinical embolization was not carried out in an optimal way for patient P. Figure 4 shows the distributions of blood concentration S 1 (t, x) and S 2 (t, x) along the AVM for the found two-stage optimal embolization regime for patient K. For patients S and P, such distributions have a similar behavior.

Discussion
Despite the well-developed technique of embolization operations, the risk of perioperative vascular rupture is still a serious danger. In this regard, the development of new methods of mathematical modeling of the AVM embolization process remains an urgent task. For modeling, it is desirable to have detailed information about the AVM geometrical structure, which would allow you to define the functions A(x) and m(x). However, modern medical examination methods of reconstructing the AVM geometry, such as computed tomography, magnetic resonance imaging and cerebral angiography, allow us to determine in vivo vessels with an average diameter of at least 0.5 mm. This not enough to determine the detailed geometric structure of the AVM, which often consisting of a large number of conjoined intertwined thin vessels, the diameter of which can reach 0.1 mm. In addition, important information is provided by intraoperative intravascular measurements of hemodynamic parameters near the pathology. Such measurements are possible only in sufficiently large vessels adjacent to the AVM, and measurements directly inside the AVM are seriously difficult and often impossible. Due to the above limitations, the proposed mathematical model is simplified, but despite this, it demonstrates a good agreement between the behavior of the calculated pressure and the clinical one. It should be noted that the question of the theoretical justification of the existence and uniqueness of the optimal control problem solution under consideration requires additional study and can be the subject of another mathematical investigation. Since this theoretical study is not the subject of interest of our applied work, we obtain optimal control using numerical modeling. www.nature.com/scientificreports/ On the basis of the mathematical model, it was possible to obtain a good correlation of the calculated and clinical pressure at the AVM input for patients K and S. Among the considered patients, patient P stands out, for which three-stage embolization was required to achieve the functional value Ĵ N < 0.05 , in contrast to patients K and S, for which such functional values were achieved in two stages. It should be noted that the patient P had a hemorrhage immediately after the neurosurgical operation, while the other patients did not have perioperative complications. This suggests that the properties of AVM and its surrounding vessels for patient P are significantly different from those of the other two patients. This is also evidenced by a significant difference in the behavior of the clinical and calculated pressure for this patient. It can be hypothesized that the need for a large number of embolization stages indicates a high probability of perioperative complications, and in this case, it is especially useful to find the optimal embolization regime. The appearance of optimal discontinuous controls ("bang-bang" control type) at the last stage for patients S and P is also interesting.
In the considered model, the influence of surrounding healthy vessels on blood flow through the AVM is taken into account by setting the q b and p 2 functions based on intraoperative measurements near the pathology. Further development of the proposed AVM embolization model may be associated with a more detailed description of the interaction of a healthy vascular network and AVM. This, in particular, will allow us to study the changes in intracerebral blood flow during embolization away from the pathology. It is also important to study the AVM, which consists of several compartments with different properties. To do this, it is necessary to improve the proposed model of optimal AVM embolization in order to take into account more physiological features of the circulatory system. Another important aspect in improving the model may be taking into account the increase in the embolic agent viscosity when it interacts with the blood.

Conclusion
In this paper it is formulated and investigated the optimal control problem of multi-stage embolization, which arises in the simulation of neurosurgery on pathological blood vessels of the brain. This simulation is performed by describing the joint filtration flow of two phases: blood and embolic agent that thromboses the pathology, as well as taking into account the redistribution of blood into healthy vessels. As a result, an initial-boundary value problem arises for a system of integro-differential equations that includes a hyperbolic divergent partial differential equation with a nonconvex flow function. The numerical solution of this problem, obtained using a www.nature.com/scientificreports/ monotone modification of the CABARET scheme, is used to solve the optimal control problem of multi-stage embolization by a modified particle swarm method.
In the numerical solution of the considered problem, the parameters and functions included in the system of integro-differential equations and the phase constraints of the optimal control problem are determined on the basis of data obtained for real patients during neurosurgical embolization operations. For the obtained solutions to the optimal control problem of multi-stage embolization, a good agreement between the calculated and clinical pressures is observed. The proposed mathematical model and numerical algorithm will be used to improve the methodology and the safety of neurosurgical operations.