{'fraction_non_alphanumeric': 0.04359834286841586, 'fraction_numerical': 0.01624251989215493, 'mean_word_length': 4.494219653179191, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Positivity preserving high order schemes for angiogenesis models

Hypoxy induced angiogenesis processes can be described coupling an integrodifferential kinetic equation of Fokker-Planck type with a diffusion equation for the angiogenic factor. We propose high order positivity preserving schemes to approximate the marginal tip density by combining an asymptotic reduction with weighted essentially non oscillatory and strong stability preserving time discretization. We show that soliton-like solutions representing blood vessel formation and spread towards hypoxic regions are captured.


Introduction
Angiogenesis is a physiological process through which new blood vessels originate from pre-existing ones [8,19], through sprouting and splitting. Blood vessel development is a normal and vital process in tissue growth, development and repair. However, it is also essential in the transition of tumors from a benign to a malignant state [42] and in many inflammatory and immune diseases. Some treatments rely on the use of angiogenesis inhibitors [9,19,33,47]. Angiogenesis is related to hypoxia, that is, low levels of oxygen. Cells in hypoxic regions release sustances, such as VEGF (Vascular endothelial growth factor), which diffuse reaching adjacent blood vessels. VEGF has been shown to increase the number of capillaries in a given network. In its presence, endothelial cells proliferate and migrate, eventually creating tubular structures resembling capillaries [30]. Fig. 1 illustrates the process.
Many models dealing with different features of angiogenesis processes have been proposed, see [24,28,29,34,41] for instance. Nevertheless, new experimental observations inspire novel mathematical models focusing on new features. The ultimate goal is to control the angiogenesis process to help healing as well as healthy tissue development. Kinetic models allow us to gain insight on the dynamics of blood vessel networks by means of asymptotic and numerical tools. We consider a kinetic model which describes the stochastic nature of blood vessel branching, derived from stochastic representations of the formation of blood vessel networks by means of ensemble averages [44]. Let us denote by p and C the density of actively moving blood vessel tips and the concentration of angiogenic factor released by hypoxic cells, respectively. Their time evolution is governed by the following system of nondimensional equations: where The term α(C)δ σv p represents the creation of new tips, while −Γp t 0 ρds quantifies tip destruction when existing vessels merge. The Fokker-Planck operator represents vessel branching and spread. The presence of the force term F(C) in the equation (1) for the tip density p forces spread towards the region where the source of growth factor is located. Similarly, the nonlocal sink −C j in (2) represents consumption of growth factor by the created vessel tips. The dimensionless parameters β, Γ, κ, χ, A, Γ 1 , δ 1 , η, and q 1 are positive. δ σv is a regularization of a delta function, a Gaussian centered about a reference velocity v 0 : An existence, stability and regularity theory in bounded and unbounded domains is developed in [11,10]. Approximate soliton-like wave solutions have been found for the density of tips, which emerge at a pre-existing blood vessel and advance until the hypoxic region is reached [4,5]: where c denotes the propagation speed and X(t) the position of the maximum. Under some hypotheses, the parameters K and c, satisfy a complementary set of equations. We focus here in two dimensional settings, where these models may be adapted to describe angiogenesis disorders causing retinopathies [7], for instance. In three dimensions, the nonlocal terms must be modified to account for three dimensional peculiarities of blood vessel interaction and high order schemes would be necessary to reproduce such complex interactions. This task may involve a large computational cost due to the presence of both space and velocity variables. Since the solutions of the model represent densities and concentrations, it is essential to use numerical discretizations ensuring positivity, a general issue with kinetic descriptions, which is challenging due to the presence of nonlocal terms and sources whose sign may change [5].
Finite difference and finite element schemes for Vlasov-Fokker-Planck models are considered in [1,2,25]. Note that system (1)-(5) resembles Vlasov-Poisson-Fokker-Planck (VPFP) equations except for the coupling with a diffusion problem and the presence of a nonlocal in time sink. In higher dimensions, random particle methods [22] and deterministic particle methods combined with splitting of diffusion and convection operators [23], have been proposed, as well as finite difference schemes combined with a change of variables [46]. These schemes are usually tested tracking convergence to a stationary solution, under periodic boundary conditions to simplify.
High order schemes for multidimensional VP and VPFP problems include approaches based on finite volume methods [3], semi-lagrangian discretization [37] or meshfree techniques [16]. In the slab geometry we consider for (1)-(5) finite differences constitute a robust choice to preserve positivity and tackle the nonlocal in time coefficients, as well as the coupling with the diffusion equation. Our goal here is to devise high order schemes. To lower the computational cost, we consider the limit of large β and work with an asymptotic approximation of the full model which provides reduced equations for key magnitudes: the marginal tip density and the angiogenic factor concentration. We will focus on the final system, which does not depend on the velocity, combining positivity preserving weighted essentially nonoscillatory (WENO) discretizations of the differential operators in space and adequate strong stability preserving (SSP) discretizations in time [48,38,39,32].
The paper is organized as follows. Section 2 presents the two dimensional reduction of the full model. Sections 3 and 4 detail the spatial discretization by means of WENO and SSP schemes. Section 5 displays numerical solutions representing blood vessel tips migrating towards a hypoxic region. Finally, Section 6 contains our conclusions.

Asymptotic reduction
in the limit of large friction, β → ∞, the kinetic equation can be approximated by a reduced equation for the marginal density (5). The source terms in equation (1) (two first terms on its right-hand side) favor velocities in a small neighborhood of v 0 . Such velocities are those for which the tip creation term proportional to α(C)δ σv (v − v 0 ) may balance the nonlocal tip destruction sink term −Γp t 0 ρds. A reduced equation for ρ is obtained by the Chapman-Enskog method [4].
In the limit as β → ∞, the marginal density ρ(x, t) and the concentration C(x, t) obey the equations [5]: To leading order, the marginal density and the density are related by A positivity preserving order one scheme would follow by explicit forward time discretization, upwind treatment of transport terms, and centered schemes for Laplacians. Integral terms can be discretized using composite Simpson rules.
To obtain a higher order scheme, we resort to positivity preserving WENO5 schemes for spatial operators, introduced in [48] for conservation laws with nonlinear fluxes and divergence free advection. Reference [48] starts from Lax-Friedrichs fluxes, we will use upwind fluxes here instead, and adapt the schemes to advection terms not necessarily divergence free. Also, the integral I(x, t) = t 0 ρ(x, s) ds becomes an additional equation The final system is composed of two equations involving diffusion, reaction, and convection terms, coupled to this additional ordinary differential equation. In the next section, we detail the spatial discretization procedure.

Weighted essentially non oscillatory space discretization
Let us consider a two dimensional convection-diffusion problem of the form where d is a constant, either d > 0 or d = 0, and c = (a, b). In our case, c = F = (F 1 , F 2 ) depends on x = (x, y) and it is not necessarily divergence free. Let us introduce a uniform rectangular mesh (x i , y j ), with x i+ 1 2 = x 1 2 + i δx, i = 0, . . . N x , and y j+ 1 2 = y 1 2 +j δy, j = 0, . . . N y . Let us denote the midpoints by x i = . The spatial steps are given by δx = X Nx , and δy = Y Ny for intervals [0, X], [0, Y ]. Several schemes for one dimensional conservations laws, and for two dimensional nonlinear laws as well as advection equations with divergence free non smooth velocities are described in [48]. We specify here the procedure we use for our two dimensional operators with non necessarily divergence free advection coefficients, adapted from [48] with upwind fluxes, instead of Lax-Friedrichs.
Step 1. Integrate (14) in the cell To handle the possible presence of a diffusive term, double cell averages are used: Integrating we find: Step 2: Replace the integrals by quadrature rules of 5-th order. Three-point Legendre-Gauss quadrature rules are the standard choice. The weights and nodes are: with γ = 1, 2, 3. For these double integrals where the final nodesx γ i and weightsw γ are with γ = 1, . . . , 5. Using these quadrature rules, (16) becomes Step 3. Upwinding in the advection terms. To ensure stability and preserve positivity, it is convenient to discretize the advection terms in such a way that information follows the characteristics [48,5]. That is achieved rewriting the terms au(x β ) in terms of numerical fluxes that are consistent with the physical flux. In the absence of diffusion and sources, this choice would also ensures that the scheme is conservative: speeds and possible shocks would be correctly captured. In this paper we use generalize upwind fluxes. Setting f (u) = au, we replace au(x β , ·) in cell C ij . Alternatively, we may use the Lax Friedrichs flux: We denote byĝ the equivalent fluxes when a is replaced by b and the roles of x and y are interchanged. In this way we obtain the scheme Step 4. Integrate in time. The time discretization strategy for (23) defines the order of the scheme in time. Forward Euler schemes are a simple choice leading to order one, which provides stability and positivity. Known the right hand r ij side of (23) at time t n , the double average u ij (t n+1 ) = u ij (t n ) + δt r ij (t n ). The next steps allow to approximate the values of u at the nodes at time t n+1 and the process is repeated at subsequent times. Higher order time discretizations must belong to the SSP class (strong stability preserving) as discussed in the next Section.
Step 5. The reconstruction of the required approximations of u from the double averages u ij , is done dimension by dimension, following a procedure detailed in [48]. We recall the procedure in an Appendix for ease of the reader. Additionally, slope limiters can are enforced as described in [48], if necessary. The scheme is expected to be stable when

Strong stability preserving time discretization
To maintain positivity and stability, we work with strong stability preserving (SSP) time discretizations. High order time discretizations of Runge-Kutta type increase the cost since several evaluations of the right hand side are required per time step. Multistep schemes avoid this problem reusing previous step, but increase storage needs and usually require smaller time steps to ensure stability.
In spite of their theoretical order, WENO5 schemes might degenerate to order two in practice in smooth regions. Therefore, we will use SSP schemes or order 2 or 3. Usual choices for third order accuracy are a third order SSP multistep method [48] u(t n+1 ) = 16 27 (u(t n ) + 3δt r(u(t n ))) + 11 27 u(t n−3 ) + 12 11 δt r(u(t n−3 )) , (26) and a third order Runge Kutta method [39] u (1) = u(t n ) + δt r(u(t n )), The stability of SSP discretizations is governed by a CFL number c as follows. If Euler forward time discretization applied to equation u t = r(u) remains stable under the condition δt ≤ δt 0 , then SSP time discretization is stable when δt ≤ d δt 0 . For the multistep method d = 1/3 while d = 1 (d ef f = 1/3) for the RK3 (27). For second order approximations, the RK2 scheme is u (1) = u(t n ) + δt r(u(t n )), u(t n+1 ) = 1 2 u(t n ) + 1 2 u (1) + 1 2 δt r(u (1) ), with d = 1 (d ef f = 1/2). In our case, the spatial operator r(u(t n )) is the operator obtained discretizing the space variables, time excluded. The final scheme for our system would read as follows.

Numerical tests in two dimensions
In this section, we present numerical solutions of the reduced equations (8)- (13) for appropriate values of the parameters as listed in Table 1. The slab geometry of the advancing angiogenic network is schematized in Figure 1. More precisely, we consider a slab (0, 1)×R, and set x = (x 1 , x 2 ). To ease the implementation of WENO discretizations, we slightly modify the boundary conditions proposed in the original model [5]. For numerical purposes we truncate the spatial domain and set Ω x = [0, 1] × [−1.5, 1.5]. We impose zero Dirichlet boundary conditions for C at x 2 = ±1.5 and x 1 = 0, as well as zero boundary conditions for ρ. At x 1 = 1, we use the function to generate boundary values. Here, c L (t) > 0 to represents the influx of angiogenic factor produced by the core of the hypoxic region and 1/a is a characteristic length thereof. This function will decrease with time as new blood vessels reach the tumor. We set in our simulations a = 1/0.3 and c L (t) = 1.1. Fot the initial data, we consider  Figures 2 and 3 show the evolution of the marginal tip density ρ and the angiogenic factor concentration. Note that we only plot the density of tips, not the full vessels already formed that lie behind them, between the pre-existing vessel and the new tips.

Conclusions
We have introduced a high order positivity preserving scheme for asymptotic reductions of kinetic models of angiogenesis, combining positivity preserving high order WENO discretizations in space and SSP discretization in time. The resulting schemes reproduce traveling patterns representing the formation and spread of blood vessel tips in the direction of an hypoxic region, in need of oxygen supply. This approach provides building blocks to devise high order positivity preserving schemes for the full kinetic models, by resorting to fractional step splitting schemes [26,18]. The proper implementation of such schemes is at present a computational challenge, out of the scope of the present study. 7 Appendix: Reconstruction from averages for ε > 0 small (typically ε = 10 −6 ) and m = 0, 1, 2. If any of the weights is negative at x, we define ω m (x) as follows for θ = 3 and m = 0, 1, 2. Once this is done, the polynomial approximation of u in S obtained from the double averages is The procedure to reconstruct a two dimensional function at the required points is as follows: • For u(x β i+ 1 2 ,ỹ α j ) − and u(x β i− 1 2 ,ỹ α j ) + -Fixỹ α j , α = 1, . . . , 5, and apply (32)-(37) to reconstruct one dimensional double averages u i (ỹ α j ) in the x direction from the values u i,j−2 , u i,j−1 , u i,j , u i,j+1 , u i,j+2 .
• For u(x α i , y j )