Blood Clotting Decreases Pulmonary Circulation during the Coronavirus Disease

: Spontaneous blood clotting in pulmonary circulation caused by thrombo-inﬂammation is one of the main mortality causes during the COVID-19 disease. Blood clotting leads to reduced pulmonary circulation and blood oxygenation. Lung inﬂammation can be evaluated with noninva-sive diagnostic techniques. However, the correlation of the severity of the inﬂammation with the pulmonary blood ﬂow has not been established. To address this question, in this work, we develop a multiscale model taking into account the interaction of a local model of thrombus growth with 1D hemodynamics in a vessel network. Flux reduction depending on the level of lung obstruction is evaluated. In particular, the model obtains that an obstruction level of 5% leads to a 12% reduction of blood ﬂux. The suggested approach can be used to investigate the interaction of blood clotting and ﬂow not only in the pulmonary network but also in other complex vessel networks.


Introduction
The coronavirus disease provoked by the SARS-CoV-2 infection is characterized by an increased risk of spontaneous blood clotting leading to deep vein thrombosis [1], cerebral thrombosis [2] and disseminated clotting in the lungs [3]. COVID-19 patients show elevated concentrations of the coagulation factors, and reinforced platelet aggregation and adhesion [4]. Various inflammatory cytokines and neutrophil extracellular traps (NET) observed during the coronavirus disease upregulate blood clotting [5]. This is particularly dangerous in pulmonary blood circulation where obstructive thrombosis can lead to a reduction in the oxygenation level with some severe cases where even mechanical lung ventilating can appear inefficient. Thus, evaluation methods of the degree of thromboinflammation and the obstruction of pulmonary circulation during the COVID-19 disease are required. Some of these methods can possibly rely on chest computer tomography (CT), allowing the diagnosis of various disorders, including obstructive thrombosis [6,7]. Lung inflammation can lead to endothelial dysfunction, under which endothelial tissue begins to express the tissue factor and initiates blood coagulation. Blood vessel obstruction from growing clots influences the pressure distribution and blood flow velocity in pulmonary circulation. On the other hand, the flow velocity regulates clot growth, and its decrease can reinforce clot formation. Altogether, inflammation and clotting lead to decreased pulmonary blood circulation and blood oxygenation. The connection between obstructive thrombosis in the lungs and the blood flow in pulmonary circulation and thus the level of oxygenation are most likely impossible to investigate with the clinical methods on the real patients.
Thus, mathematical modelling can be useful here. However, 2D and 3D modeling of blood flow in a complex network of vessels cannot be used because of computational limitations. Therefore, the quasi-one-dimensional approach [8,9] is widely used in this case.
In this work, we develop a novel mathematical model describing the interaction of blood clotting in small lung blood vessels with pulmonary circulation. The quasione-dimensional model is used to describe pulmonary circulation, and a single reactiondiffusion equation models blood coagulation. The 2D flow model allows us to determine the flow velocity at the clot surface from the results of quasi-one-dimensional simulations. The quasi-1D model is derived under the assumption that the vessel length is much larger than its diameter and that the flow is radially symmetric. The latter assumption is satisfied everywhere except the clot area. A comparison with the direct numerical simulations of the 2D problem shows that the 1D approximation is applicable in the clot area also.

Hemostasis and Thrombo-Inflammation
Blood coagulation is an important physiological mechanism that aims to stop bleeding in the case of injury. It includes fibrin gel formation in the blood plasma and platelet aggregation. These two processes interact with each other, and they are regulated by endothelial cells and blood flow.
In the case of injury, the tissue factor (TF) in the subendothelial layer of the blood vessel wall comes into contact with blood plasma and initiates the coagulation cascade. Factor VII in blood plasma forms a complex with TF and activates factors IX and X, initiating the conversion of prothrombin into thrombin. If the initial thrombin concentration produced near the vessel wall exceeds a certain threshold level, it starts the self-amplifying loop of thrombin production by means of factors V, VIII and XI ( Figure 1). Thrombin is an enzyme that converts fibrinogen in blood plasma into fibrin, which then polymerizes and forms a fibrin clot. On the other hand, the tissue factor pathway inhibitor (TFPI), antithrombin and activated protein C act as anti-coagulant factors, whose roles are to control or stop the coagulation cascade. In the case of injury, the tissue factor initiates thrombin production. Thrombin is an enzyme that converts fibrinogen into fibrin, leading to the formation of fibrin clot. Thrombin (IIa) catalyzes the activation of factors V, VIII and XI; factors XIa and IXa catalyze the activation of factors IX and X, respectively; factors VIIIa and Va form active complexes with factors IXa and Xa, respectively, and further increase thrombin production.
Platelets are small blood cells whose concentration increases near the vessel walls due to margination by red blood cells. Being activated by thrombin, some other molecules, other activated platelets or a high shear rate, they begin to aggregate and express various pro-coagulant factors, upregulating the coagulation cascade in the blood plasma. Thus, coagulation reactions and platelet aggregation work together to stop bleeding. Platelets' aggregation prevails in the case of arterial clots (white clot), while fibrin generation plays the most important role in veins (red clot).
Blood coagulation represents a subtle equilibrium between pro-and anti-coagulant factors. The violation of this equilibrium can lead to insufficient coagulation and bleeding disorders (e.g., hemophilia), while excessive coagulation results in thrombosis, possibly leading to thromboembolism, heart attacks or ischemic strokes.
Hemostasis can be influenced by various pathologies and physiological conditions, such as cancer, diabetes or inflammation [10]. As an example, rheumatoid arthritis leads to elevated concentrations of TNF-α and IL-17 in blood, resulting in a 10-fold increase in the risk of cardiovascular events related to thrombosis and atherosclerosis [11].
SARS-CoV-2 infection influences blood coagulation in multiple ways by upregulating the coagulation cascade and platelet aggregation. Pro-inflammatory cytokines produced during infection development, such as TNF-α, IL-1β, IL-17 and many others, promote endothelial dysfunction, initiating the blood coagulation cascade at the vessel walls [5]; TNF-α inhibits thrombomodulin expressed by the endothelial cells whose complex with thrombin activates protein C, which provides one of the main mechanisms of clot growth arrest; IL-17 inhibits CD39, which is an inhibitor of platelet aggregation [11]. Platelets in COVID-19 patients show increased activation and aggregation. Immune cells, neutrophils, produce various pro-inflammatory cytokines (likewise macrophages), reactive oxygen species (ROS) and neutrophil extracellular traps (NET), all of which promote thrombosis. Thus, thrombo-inflammation in COVID-19 represents one of the main risk factors.

Combining Local Models of Blood Clotting with Global Models of Blood Circulation
Blood coagulation, initiated by damage to a blood vessel wall, leads to the formation of a fibrin polymer and propagates inside the vessel as a reaction-diffusion wave. Clot growth is influenced by blood flow and, in turn, changes the flow velocity and pressure distributions. There are numerous modeling studies of clot growth in quiescent plasma and in flow (see [12] and the references therein). All these works are devoted to the investigation of a small part of a blood vessel in the vicinity of a growing clot. The question about its interaction with global blood circulation was not previously studied.
The aim of the present work is to develop a model of multiple clot growths in the vascular network while taking into account their mutual influence on blood flow in the network. Such model includes three parts: a thrombus growth model, a blood flow model and a model of the vessel network. The thrombus formation model and the model of blood flow are interrelated: a growing thrombus changes the flow conditions, while blood flow carries away blood factors and prevents thrombus formation. To study this process in pulmonary blood circulation, a model of such network is required.
The reduction to the 1D problem of clot growth is done without the symmetry assumption in the following way. We consider the clot in the vessel in the 2D formulation [5], and there is no symmetry assumption. With this model, we obtain the speed of the thrombin propagation wave [13]. This speed determines whether the clot is growing. If this is the case, the clot reduces the cross-section area of the vessel. This reduced cross-section area of the vessel is then used in the quasi-one-dimensional system. The results given by this approach were tested on the network consisting of the three vessels (Section 3.4), and the results agree with the 2D calculations of the same problem.
In the next section, we present a simplified model of clot growth in terms of a reaction-diffusion equation for thrombin concentration. This model, validated in previous works [14], allows us to determine the final clot size depending on the flow velocity. Section 3 is devoted to models of blood flow. We use a 1D hemodynamic model to describe the flow velocity and pressure distribution in a vessel network. This model is combined with 2D Poiseuille flow in order to restore the flow velocity distribution in the vicinity of the clot. This approach is validated by a comparison with 2D numerical simulations of the Navier-Stokes equations. The model of pulmonary vessel networks is presented in Section 4. Finally, in Section 5, we bring all these models together and describe blood clotting in a vessel network. We determine the degree of vessel occlusion and the flux through the network depending on pressure difference and the rate of blood coagulation.

The Model of Clot Growth
Blood coagulation in a quiescent plasma can be modeled with reaction-diffusion systems of equations for the concentrations of blood factors [13][14][15]. Clot growth, that is, the spatial progression of the region filled by fibrin polymer gel, corresponds to the propagation of a reaction-diffusion wave [16]. A simplified model of blood coagulation consisting of three equations for the concentrations of prothrombin, thrombin and activated factor X was developed in [17] with patient-specific parameters. These results are used in the present work to determine the parameters of the one-equation model for the timedependent thrombin concentration T: Here, y is the spatial coordinate in the direction perpendicular to the vessel wall, t is time, and D is the diffusion coefficient. Function Φ(T) represents the thrombin production rate. It takes the following form after the kinetic model reduction [14,16]: where A 0 is the antithrombin concentration assumed to be constant, P 0 is the prothrombin concentration initially present in the flow, k 3 , k 4 , k 5 and k 6 are kinetic constants, α is a parameter specified below, and u is the blood flow velocity. Let us note that the convective term u ∂T ∂x in the axial direction can be approximated by u (T − T 0 )/L, where T 0 = 0 is the thrombin concentration at the entrance of the clotting area, and L is its length. This gives us the term −αuT.
Generally speaking, the blood flow velocity u is a function of the spatial coordinate in the vertical direction. However, it can be considered constant at the clot surface. Previously, this model was considered in [14,16] in order to determine the speed of clot growth and the final clot size. This simplified model gives a good approximation of more complete models, including the model of clot growth in flow with Navier-Stokes equations for the blood flow velocity.
Equation (2) considered on the whole real axis has a traveling wave solution T(x − ct), satisfying the problem Here, c is the wave speed corresponding to the speed of clot growth, and T * is the maximal zero of the function Φ(T). There are different methods to determine the speed c, e.g., [16]. However, here, we are interested in the condition under which the clot stops growth, that is, c = 0. This condition can be derived by the following analytical method. We multiply Equation (3) by T and integrate from minus to plus infinity: where ξ = x − ct. The first integral with the boundary conditions from Equation (3) equals zero [18]. From the last two terms, we obtain The denominator of this expression is positive, so the sign of the wave speed is determined by the integral in the numerator. If this integral is positive, then the clot grows; if it equals zero, then the clot stops growing. We will use this condition below to determine the final clot size.
Let us now estimate the blood flow velocity u for which clot growth stops. We use the equality and the values of parameters adjusted to [17] corresponding to the normal level of blood coagulation ( Table 1). Calculations of the one-equation model with these parameters without blood flow gives the wave speed of thrombin propagation 0.54 micron/sec corresponding to the physiological range.  (1)-(2) for the normal blood coagulation estimated from [17]. The speed of thrombin concentration wave propagation for these parameters is 0.54 micron/s, which corresponds to the evaluative value for the normal blood coagulation [17]. The speed is calculated under the assumption of the absence of blood flow (u = 0).

Parameter Value
Diffusion coefficient D 5×10 −5 cm 2 /min Initial concentration of prothrombin P 0 1400 nM Under the approximation T * ≈ P 0 , we have where T = θP 0 , σ = k 6 A 0 + αu and the terms with k 3 and k 4 are neglected since they are small in comparison to the next term. Equating the last integral to 0, we obtain In order to simplify the interpretation of this relation, let us suppose that k 6 A 0 αu. This assumption is satisfied for the realistic values of parameters. Taking into account that α ∼ 1/L, where L is the clot length, we conclude that αu = 1/τ f , where τ f is the characteristic time that blood flow stays in the clotting area. Next, k 5 P 3 0 = 1/τ r , where τ r is the characteristic reaction time. Therefore, we obtain the relation between the reaction time and residence time that provides the clot growth arrest. If the blood flow velocity is small, such that the residence time exceeds the required threshold, the clot continues to grow.
Let us note that the flow velocity u is considered here as constant. In what follows, we will consider the Poiseuille flow where u = u(y) is a parabolic function. In this case, instead of the exact condition in Equation (5), we will use an approximate condition considered for a fixed value of y, where This approximation is justified since the reaction zone of the thrombin wave, that is, the interval where thrombin concentration essentially changes, is sufficiently narrow, and the variation in the flow velocity there is negligible. Figure 2 shows the solution of Equation (1) with the space-dependent nonlinearity Φ(T, y) due to the parabolic flow velocity profile u(y). The thrombin wave propagates until the flow velocity becomes sufficiently large enough to stop it. If the maximal flow velocity at the center of the vessel is not sufficient to stop the wave, then it propagates till the end of the interval (Figure 2, right). This case corresponds to complete vessel occlusion.  Table 1, α = 0.5 for the left image and α = 0.4 for the right image, H = 1 cm.
For the parameter values in Table 1 for normal blood coagulation and the clot length 1 cm, we obtain the critical flow velocity of the order of 0.1 cm/min from Equation (6). If the flow velocity is larger than this value, the clot does not grow. We will see below that the average flow velocity in vessels of generations 9-12 of pulmonary circulations is about 10 3 times larger. However, the flow velocity near the vessel wall (where the clot starts growing) is essentially less. Moreover, the value of the constant k 5 increases in the case of COVID-19 patients. Finally, the inflamed area can also be sufficiently large. These factors taken together can lead to the emergence of inflammation-induced clotting in pulmonary circulation and to the decrease of the total blood flow.

The Model of Blood Flow
The goal of this work is to investigate the mutual influence of blood flow and clot formation in a network of pulmonary blood vessels. Due to the complexity of this network, we will consider simplified flow models with Poiseuille flow and 1D hemodynamics. The flow velocity distributions will be compared with the solution of the Navier-Stokes equations for the model problem in three connected vessels.

Navier-Stokes Equations
We begin with a single vessel with the clot schematically shown in Figure 3 [14]. In order to compare numerical results with the analytical solution considered below, we approximate clot shape by the rectangular domain. The flow in this vessel is described by the Navier-Stokes equations: Here u is the flow velocity, p is dynamic pressure, ρ is density of the blood, µ is the dynamic viscosity, and K f is the hydraulic permeability of the clot [19]: where α is the fiber radius. Thus, we consider the whole rectangular domain, including the clot that modifies the flow velocity by the additional volume force. No-slip boundary conditions for the velocity are imposed on the upper and lower boundaries, and parabolic velocity profiles at the entrance and exit of the vessel. These parabolic boundary conditions are imposed up to a constant factor that is determined from the condition of a given pressure difference between the entrance and the exit of the vessel. Such boundary conditions allow complete vessel occlusion by growing a clot contrary to the conventional Dirichlet boundary condition with a given flow velocity at the boundary of the domain [14]. Blood flow is assumed to be Newtonian, and we perform numerical simulations using the Computational Fluid Dynamics library OpenFOAM [20]. We will compare below the result of numerical simulations of this model with the analytical approximation in the case of a single vessel presented here. Then, we study the case of bifurcating vessels.

Poiseuille Flow
Consider the flow in the rectangular non-deformable vessel with length L and height H without a clot. Suppose that the transverse velocity component u y equals zero, and longitudinal velocity component u x is parabolic as in the Poiseuille flow: The parameter a is unknown. In this case, we obtain from the Navier-Stokes equations: Using the assumption of the parabolic velocity profile from Equation (10), we obtain the following expression for the parameter a: where ∆p is the pressure difference between the entrance and the exit of the vessel. The average u av and the maximum u max velocities are given by the following expressions: Let us now consider the flow in a vessel with a rectangular clot with width l and height h (Figure 3, right). We approximate the flow velocity by parabolic profiles, which are different above the clot and outside the clot. The velocity above the clot vanish not at the clot surface y = h but below it at y = h − d: This assumption implies nonzero flow velocity at the clot surface, which corresponds to the experimental observations and allows us to describe the influence of blood flow on clot growth. From Equation (10) and mass conservation, we obtain: From the condition that the total pressure difference equals ∆p, we find: where k is given in Equation (15). Average and maximal velocities before and after the clot are determined by Equation (13), where a is defined by Equation (16), while the average velocity above the clot and the velocity at the clot surface can be determined by the formulas: where a is given by Equation (16) and k by Equation (15).

Blood Flow in the Network of Three Vessels
Let us now consider the network of three vessels (Figure 4). We suppose that there is a rectangular clot with the height h and the length l in one of the branching vessels. As before, we assume that the flow velocity has the parabolic profile in the parts of the network without a clot, and it has a "shifted" profile as in Equation (14) above the clot. The mass and momentum conservation are imposed at the bifurcation point. Note that a similar configuration of the three vessels with a clot was considered in a recent work [21], but the clot was considered as completely impenetrable for fluid flow. This difference is important for the present study because the flow velocity at the clot surface determines the speed of its growth. Using Formulas (12) and (16), from the mass conservation, we obtain: where l i and h i are the length and the height of the vessel i, respectively, i = 0, 1, 2, p 0 and p 2 are the left and the right boundary pressure, respectively, p 1 is the pressure at the bifurcation point, l is the clot length, and parameter k is defined by Equation (15). Solving this equation, we obtain the pressure p 1 at the bifurcation point: We compare the analytical solution with numerical simulations of the Navier-Stokes Equation (8) (Figure 4, right). We suppose that the permeability of the clot corresponds to a fiber radius equal to 0.2 µm and a uniform concentration of the fibrin polymer of 4900 nM. The difference in the pressures between the inlet and the outlets is kept constant and equals 30 g/ cm × s 2 . Other parameters used in the simulation are listed in Table 2.
The analytical solution (black lines) shows good correspondence with the numerical results (red lines) ( Figure 5A-C). The influence of clot growth on hemodynamics is studied by varying the clot height and evaluating the pressure ( Figure 5A) and velocity ( Figure 5B) at the center of each vessel. The results show that the velocity at the vessel 0 (before the bifurcation point, Figure 4) slightly decreases as the clot height increases, while the pressure decreases.
The effect of clot growth on hemodynamics is different in peripheral vessels 1 and 2. Clot growth decreases the flow velocity and increases the pressure in vessel 2, where the clot is located ( Figure 5). However, both the velocity and pressure increase in vessel 1. The pressure at the bifurcation point also increases with the increase of the clot height.
To evaluate the effect of clot growth on flow velocity above the clot, we increase the height of the clot, and we measure the maximal value reached by velocity in the space between the clot and the vessel wall. Numerical and analytical results are shown in Figure 5C. In both cases, the maximal value of the velocity changes depending on the thrombus size. This maximal value increases as the clot height increases. As the clot height reaches some value (550 µm in numerical simulation), the velocity value starts decreasing until it reaches zero when the clot completely occludes the vessel. Pressure p 2 at the right end of the network 0 g cm·s 2 Coefficient d from Equation (14) 0.003 cm The analytical approach presented here will be used below in more complex networks of blood vessels in order to determine the flow velocity at the clot surface from the average flow velocity obtained in the model of 1D hemodynamics.

Quasi-One-Dimensional Flow
Quasi-one-dimensional approximation is widely used in modeling the blood flow in complex systems of vessels [8,[22][23][24]. The model consists of the following equations [23,24]: where u(x, t) is blood velocity, p(x, t) is pressure, s(x, t) is a cross-section area of the vessel, ν is kinematic viscosity, ρ is density, x is spatial coordinate, x ∈ (0, l), l is the vessel length, and t is time. This model represents mass and momentum conservation equations. The right-hand part of the momentum equation corresponds to the Poiseuille flow if the cross-section area is constant. In the case of several vessels, mass and momentum conservation provide the boundary conditions at the bifurcation points. This system is completed with the third equation describing the pressure-area dependence. For the vessels without the clot, we will use the following pressure-area dependence [23]: where s 0 , θ and p 0 are parameters. Some examples of other pressure-area relationships and their impact on the flow are presented in [25]. Some of the vessels of the network will be considered as inflamed where the sub-endothelial layer can initiate clot growth along the whole length of the vessel. We will suppose that the cross-section of these vessels is independent of pressure (due to larger stiffness), but it can depend on the growing clot.
To determine the final clot height and the final cross-section area of the inflamed vessels, we solve the system of equations of 1D hemodynamics and find the average flow velocity u av for the initial cross-section area without the clot. Then, the parabolic profile above the clot is restored for this average velocity, and the flow velocity at the clot surface u clot is determined: Finally, the flow velocity at the clot surface u clot is used to determine the sign of the integral (Equation (4)) and then the sign of the thrombin wave velocity c. If it is positive, the clot grows, and it leads to decreasing the cross-section area (vessel height). In this case, we recalculate the average velocity for the smaller cross section area of this vessel. If it is zero, then the clot stops growing, and its maximal size is found. Otherwise, we proceed to the next step with a decreased cross-section. Finally, either we find the cross-section for which the flow velocity stops clot growth, or it grows till complete vessel occlusion. This approach gives a stationary solution of the problem of clot growth in a network of blood vessels. We test this model on a network of three vessels (Figure 4) and compare the result with the analytical model ( Figure 6). The pressure in the quasi-1D simulation ( Figure 6A, blue line) is close to the analytical solution (black line). Average velocities in the quasi-onedimensional case are also close to the analytical solution ( Figure 6B, left). The main difference between numerical results for quasi-1D hemodynamics and the analytical results is observed for the average velocity above the clot ( Figure 6B, right). For the intermediate values of clot height, the quasi-one-dimensional approach gives larger values for the average velocity above the clot than the analytical value. As a result, the vessel is not fully occluded for the smaller pressure difference in the quasione-dimensional case than in the analytical solution (Figure 7). This difference in the estimation of the average velocity above the clot occurs in the transition from the quasione-dimensional to the two-dimensional flow field. For the simplicity of calculation, this transition is done under the assumption of mass conservation similar to the single vessel. However, in the vessel network, this condition leads to velocity difference. We calculate the correction for the considered network of three vessels: The results for the quasi-one-dimensional model taking into account this correction are shown in Figure 6 with magenta lines. They completely coincide with the analytical solution (black lines). This correction depends on the network topology, but its calculation becomes cumbersome for complex networks. Therefore, in what follows, we will use the same calculation as for the single vessel. In this case, we obtain an upper estimation for the average velocity in comparison with the exact solution.

The Network of Pulmonary Vessels
In this section, we will study the interaction of pulmonary circulation with the model of clot growth considered above.

Physiology of Pulmonary Circulation
The pulmonary circulation serves for blood oxygenation. It starts from the right ventricle and ends in the left atrium, and as a systemic circulation system, it includes arteries, capillaries and veins. The main pulmonary artery starts from the right ventricle and branches into the left pulmonary artery and the right pulmonary artery. They go into the left and the right lungs, respectively, and further branching occurs inside the lungs. There are 15 vessel generations in the arterial part of the pulmonary circulation system [26]. After passing the arterial part, blood flows to the capillaries. They are placed in the alveolar walls, where blood oxygenates. From the capillaries, it collects in the venules connected to veins. The four biggest veins are open to the left atrium. There are 15 vessel generations in the venous part of pulmonary circulation [26].
Pulmonary circulation is a low-pressure system. In the main pulmonary artery, the systolic pressure is 15-30 mm Hg [27], with an average pressure of about 25 mm Hg [28]. The diastolic pulmonary arterial pressure is varying in the interval 4-12 mm Hg [27], with an average of about 8 mm Hg [28]. The mean pulmonary arterial pressure is about 15 mm Hg [28]. The mean pulmonary capillary pressure is about 7 mm Hg [28], and such low pressure is essential for the blood oxygenation in pulmonary capillaries. The pressure in the left atrium and the major pulmonary veins is varying from 1 mm Hg to 5 mm Hg, with an average of about 2 mm Hg [27,28].
The blood volume of the lungs is about 450 mL, and it is about 9% of the total blood volume of the entire blood circulation [28]. Approximately 70 mL of this pulmonary blood volume is inside pulmonary capillaries, and the remainder is divided equally between pulmonary arteries and the veins [28].
The total flux in pulmonary circulation is the same as in the big circulation: about 80-100 mL/s.

Pulmonary Vessel Branching
There are three models of a pulmonary vessel tree [29]: Weibel model, Strahler model and Diameter-Defined Strahler's model. In the Weibel model, the tree of pulmonary vessels is supposed to be symmetric. This assumption leads to an overestimation of the total amount of vessels. In the Strahler model, the vessel tree is not symmetric. Vessel numeration starts from the smallest arterioles. At the bifurcation points, if the diameter of the next vessel is larger than the sum of the diameters of joining vessels, it is considered that this vessel belongs to the next generation. Otherwise, it is counted in the same generation. In the Diameter-Defined Strahler's model, it is supposed that the vessel belongs to the next generation if its diameter is larger than the sum of diameters of joining vessels taking into account the standard deviation of the diameters.
The Strahler model was used by Horsfield et al. to measure the diameter, the length and the number of vessels in each generation of pulmonary arteries [30] and veins [31]. In [26], these characteristics were measured with Diameter-Defined Strahler's (DDS) approach for both arterial and vein trees, and a comparison with Horsfield's data was carried out. The total cross-section area of arteries and veins for each generation was investigated. In both sets of data, the total cross-section of arteries in each generation is larger than the total cross-section of veins [26], which contrasts with the data for the big circle of blood circulation [28], where the total cross-section area of veins exceeds the total cross-section area of arteries by almost four times. In [26], this difference is even larger than in Horsfield's data. The DSS approach takes into account more factors, and the results are more accurate [26]. In what follows, we will use vessel parameters and statistical information from [26].

The Graph of Pulmonary Vessels
We use the data from [26] in order to model a graph of pulmonary vessels and to estimate the flow in the vessel network with different degrees of vessel occlusion. We will consider four consecutive generations of vessels from generation 12 to generation 9 ( Figure 8). The parameters of the graph are given in Table 3. Blood vessel diameter and length in successive generations are estimated according to Ref. [26]. The linear dependence of the cross-section area on pressure is used following Ref. [32]. The data on the maximal and minimal cross-section areas used to calculate the elasticity coefficient q are taken from Ref. [33]. The initial pressure is considered to be equal to the capillary pressure, while the initial flux is determined from the average value of 80 mL/s of normal flux in pulmonary circulation [24]. The boundary conditions for the graph are represented by the pressure values at both ends of the network estimated from the assumption of the linear pressure distribution discussed below. Figure 8. The graph of pulmonary vessels considered in the present work. The number above the vessel indicates the generation in the vessel tree [26]. The pulmonary tree starts from generation 1 corresponding to capillaries and ends with generation 15 corresponding to the main pulmonary artery leaving the right ventricle. Inflamed vessels in the 10th generation are shown in bold. Table 3. Parameters of the graph of pulmonary vessels. The diameters and length of the vessels are taken from [26]. The cross-section area is recalculated from the diameter. For all vessels, the maximal cross-section area s max and the minimal cross-section area s 0 are, respectively, 1.4 and 0.9 of the value given in the table [33]. The minimal pressure p 0 and maximal pressure p max equal 2 and 15 mm Hg, respectively [27]. In Equation (22), parameter θ = (s max − s 0 )/(p max − p 0 ). The initial pressure is 7 mm Hg for all vessels. The initial flux is 0.64 mL/s in level 12, which corresponds to the 80 mL/s in the main pulmonary artery.

Vessel
Diameter, cm Cross-Section Area cm 2 Length cm In systemic circulation, the essential part of the pressure drop occurs in precapillary vessels [29]. The pressure distribution in pulmonary circulation is reported to be more symmetrical than in the systemic circulation [34] in the sense that the pressure in pulmonary capillaries approximately equals an average between the pulmonary arterial and venous pressures. Furthermore, the pressure gradient across the vasculature in pulmonary circulation is close to constant [29]. In the current work, we consider four consecutive generations of arterioles of the pulmonary network tree. According to data from [29] (Figure 4 in the source), the pressure decay in this network can be considered linear. Therefore, the linear distribution is used to estimate the pressure at the inlet and outlet of the network.
The pressure in the pulmonary artery varies from 8 mm Hg (average diastolic) to 25 mm Hg (average systolic), and the mean pressure is 15 mm Hg [27]. In pulmonary capillaries, the mean pressure is 7 mm Hg [27]. Consider the mean pressure in the arterial part of the pulmonary vessels. It decreases from 15 mm Hg (mean) in the pulmonary artery to 7 mm Hg in the capillaries, so the pressure difference for the arterial part is 8 mm Hg, while its length is 3.62 cm (Table 3). Hence, the linear pressure distribution allows the determination of the pressure difference in the considered part of the arterial tree: 9.58 -7.75 = 1.83 mm Hg ≈ 2 mm Hg.
Next, consider the linear pressure distribution for the venous part. The total length of the venous part is 7.24 cm (calculated with the data from Huang and Yen [26]). The pressure in the left atrium varies from 1 mm Hg to 7 mm Hg, and the mean pressure is 2 mm Hg [24]. The pressure in the big pulmonary veins is supposed to be the same as in the left atrium. Therefore, in the venous part, the pressure decreases from 7 mm Hg to 2 mm Hg.
Thus, we determine the pressure at the inlet of the network from the linear approximation. The value of the pressure at the outlet of the network is determined from the solution of the quasi-one-dimensional problem in the network in such a way that the total flux through the whole network of pulmonary vessels is 96 mL/s. This value lies in the physiological range of 80-100 mL/s estimated for pulmonary circulation.

Pulmonary Circulation in the Case of Blood Clotting
The vessel network studied in this work includes four generations of arteriols ( Figure 8) with their diameter decreasing from generation 12 to generation 9. In the model, the vessels of the same generation have the same cross-section area, length and elasticity.
We proceed with the calculation of the flow in this network. The results for the network without inflamed vessels are shown in Figure 9 (left). There are different paths indicated by different colors. Each path consists of four vessels, and they differ from each other by the values of flow velocity and total flux. The plots of the pressure, velocity, crosssection area and flux have the same color as the corresponding path in the network. The pressure decreases linearly inside each vessel, but its slope changes from one generation to another one. The cross-section areas are the same for all vessels from the same generation. The difference between the paths can be observed in the values of velocities and flux. In the 10th generation, the smallest average velocity is in the orange vessels, while the largest one is in the green vessels; however, the largest average velocity is reached in the blue vessels. The total flux in the network is 0.758 mL/s. The corresponding total flux in the whole pulmonary tree is the total flux in the network, i.e., the flux in the vessel of the 12th level, multiplied by the number of vessels on the 12th level from [26]. This flux equals 96 mL/s, and it lies in the physiological range for the flux in pulmonary circulation. The volume of the represented network is 0.0748 mL.
Next, consider the network with inflamed vessels, assuming that all vessels in the 10th generation are inflamed ( Figure 8). From the biological point of view, this assumption means that inflammation provokes endothelial dysfunctions in these vessels, leading to a possible clot growth along the whole vessel length. However, whether it will grow or not depends on flow conditions: a large flow velocity prevents the clot growth, while a small velocity can lead to clot growth. On the other hand, clot growth influences the flow velocity distribution by either decreasing or increasing it depending on the geometry of the networks, pressure difference and the parameters of blood coagulation.
We repeat the same simulations in this network where a part of the vessels is inflamed. The results are shown in Figure 9 (right). At the bottom, the final clot size is evaluated. The largest clots in the 10th generation develop in the orange vessels where the average velocity is minimal in the case without inflammation. These vessels become totally occluded, and the flow velocity in these vessels and in the vessels after them is zero. Smaller clots develop in green vessels, where the velocity in the case without inflammation is higher. It obstructs up to about 11.2% of the vessel diameter. The smallest clots are observed in the blue vessels where the velocity is maximal. They obstruct about 4% of the vessel diameter. The flow velocity increases in the partially obstructed vessels. The cross-section area only reduces for the obstructed vessels, and it does not change for the other vessels. The total flux in the network reduces to 0.658 mL/s, which is 86.8% of the flux in the previous calculation without inflammation. The volume of the network after the clot develops is 0.0708 mL, with a reduction of about 5%. Though the network volume reduction as a result of clot developing is small, the flux reduction is much more significant. We will now study how the volume and flux reductions depend on the kinetic constant k 5 , characterizing the rate of thrombin production. It can be essentially increased due to inflammation. The flux dependence for two different pressure differences at the inlet and outlet of the network is shown in Figure 10A. In both cases, the flux decreases as the coefficient k 5 increases. This dependence is not linear. For smaller k 5 , the flux decrease is slow, accelerating with the increase of this coefficient. This behavior can be explained by the interaction of the growing clot and blood flow velocity in the inflamed vessels. Increasing k 5 results in decreasing the cross-section area of inflamed vessels due to clot growth, but, in turn, the velocity in the other inflamed vessels increases, and clots grow slower or completely stop growing. As a result, the total flux reduces slowly. With a further increase of k 5 , the increased velocity in the other inflamed vessels is not sufficient anymore to stop growing clots. When the clot reaches the middle of the vessel diameter, it cannot be stopped by the blood velocity since the flow velocity decreases, and the flux decreases more rapidly. This is why the flux dependence on k 5 in Figure 10A has inflection points for k 5 = 4.1 and k 5 = 5.4. When k 5 is sufficiently large, the total flux through the vessel network becomes zero, which means that all inflamed vessels are occluded.
Let us define the network volume as the total volume of all vessels with nonzero flow velocity, taking into account vessel narrowing due to clot growth. The dependence of the network volume on the coefficient k 5 is presented in Figure 10B. As before, the volume decrease also accelerates as k 5 increases due to the threshold change of the cross-section area of the inflamed vessels as the clot grows. Figure 10C shows the dependence of the flux on the volume of the network. This dependence is almost linear, and it is the same for different pressure gradients. Thus, the flux reduction can be estimated for different levels of the network occlusion. As such, only 10% of network obstruction leads to 35% of flux reduction. Figure 10. Dependence of the total flux through the network (A) and volume network (B) on the kinetic constant k 5 for different pressure gradients. Flux-volume dependence (C) is close to linear, and it does not depend on pressure difference.

Conclusions
We develop in the present work a novel method to study the interaction of blood clotting with hemodynamics in a complex vessel network. We combine a local model of clot growth with a global model of blood circulation. The model of clot growth consists of a single equation to calculate the thrombin concentration, which gives a good approximation of more complete coagulation models [14]. The flow velocity in the vessel network is modeled with 1D hemodynamics and 2D Poiseuille profiles. The latter allows restoring the flow velocity field in the vicinity of the clot. This approach is validated by the comparison with the solution of the Navier-Stokes equations.
This approach allows us to determine the flow velocity in the network of vessels with growing clots. We determine the flow velocity and pressure distributions, as well as the final clot height, total flux and network volume for different values of kinetic constants and pressure gradient.
We obtain that even small lung obstruction can lead to an essential flux reduction. In particular, 5% of the network volume reduction due to clotting leads to 12% of the flux reduction and the corresponding drop in blood oxygenation. Let us recall that an oxygen saturation level below 90% is considered as low, and below 80% as life-threatening. The value of the flux reduction depending on the network obstruction represents an estimate from above, so the real flux decrease can be less than the model prediction.
In order to model blood clotting in complex networks of pulmonary blood vessels, we have to introduce some simplifying assumptions. In particular, we assume that the vessel walls in the clotting area are rigid and that blood is a Newtonian fluid. In this case, blood flow can be described with the Navier-Stokes equations. In the analytical part of the work, we suppose that the clot has a rectangular shape and the flow velocity has a parabolic profile. Then this model admits an analytical solution, and it shows good agreement with the 2D Navier-Stokes solution. Next, we simplify the coagulation cascade and reduce it to a single equation for the thrombin concentration. This simplification allows the determination of simple analytical conditions of clot growth arrest (zero wave speed), as validated in the previous works [16]. Furthermore, since the pulmonary arterial tree contains 15 generations with millions of blood vessels, we consider a model example with only four generations of blood vessels. More complex networks can be considered in future work using the basis of the method developed in this work.
The limitations of the model are determined by the simplifications described above. Among the other limitations, we should indicate that the blood flow velocity is pulsating, which can potentially influence the clot growth rate. Furthermore, patient-specific parameters of the coagulation cascade and of lung inflammation should be taken into account for clinical applications of the developed method. Despite all these restrictions, we expect that the novel model developed in this work will open new perspectives in the investigation of thrombo-inflammation in complex networks of blood vessels.
The limitations of the approach developed in this work are mainly related to the simplifications of blood coagulation and blood flow models. However, since they are validated by comparison with more detailed models, we can expect that the results give a realistic estimation of pulmonary circulation depending on the lung damage due to inflammation. Blood flow pulsations and platelet aggregation are not considered in this work and will be taken into account in future studies. Let us also note that this modeling approach can be used for other network configurations and for other pathologies, including disseminated blood coagulation in sepsis.