α+ is a characteristic rate constant of the Ca-troponin; *ICa*(*t*) = *I* 0 *Ca*(*exp*(−*k*1*Cat*) − *exp*(−*k*2*Cat*)) is a given inflow of Ca2+ ions to the cell normalized by с0 50. The parameters of this block of the model and their values are presented in Table S1.

#### *A.2. Left Ventricle Approximation*

LV approximation was set by the following expressions for Cartesian coordinates (*r*, *z*) through the curvilinear coordinates (γ, ψ) with parameters ε*, rin, rout, hin,* and *hout* [16,22].

$$\begin{array}{l} r = (r\_{in} + \gamma(r\_{out} - r\_{in}))(\varepsilon \cos \psi + (1 - \varepsilon)(1 - \sin \psi)), \\ z = (h\_{in} + \gamma(h\_{out} - h\_{in}))(1 - \sin \psi) + (1 - \gamma)(h\_{out} - h\_{in}). \end{array} \tag{A11}$$

#### *A.3. Hemodynamics Model [17]*

Passive pressures of the right ventricle and the atria were described as follows.

$$\begin{split} P\_{LA\_{\text{Pas}}} &= E\_{1LA} \cdot \left( e^{E\_{2LA}V\_{LA}(t)} - e^{E\_{2LA}V\_{\text{OLA}}} + \mu\_{a1} \frac{\partial V\_{LA}(t)}{\partial t} \right), \\ P\_{RA\\_Pas} &= E\_{1RA} \cdot \left( e^{E\_{2RA}V\_{RA}(t)} - e^{E\_{2RA}V\_{\text{OLA}}} + \mu\_{a1} \frac{\partial V\_{RA}(t)}{\partial t} \right), \\ P\_{RV\\_Pas} &= E\_{1RV} \cdot \left( e^{E\_{2RV}V\_{RV}(t)} - e^{E\_{2RV}V\_{\text{OLV}}} \right). \end{split} \tag{A12}$$

Here, *P\*\*\_Pas* are passive elastic parts of chamber pressures, *V\*\** are their volumes, and *V*0*\*\**, *E*1, *E*2, and μ*a*1 are constant parameters. Subindices *LA*, *RA*, and *RV* stand for the left atrium, right atrium, and right ventricle, respectively. Active pressures of the right ventricle and the atria were found from ordinary differential equations. Due to the introduction of time delay with relaxation time τ and the analogue of force–velocity equation, those described the pressure time-course more accurately than the pressure–volume dependencies with time-dependent stiffness coefficients commonly used in other lumped parameter models.

$$\begin{aligned} \tau \frac{\frac{\partial P\_{LA\\_Act}(t)}{\partial t} + P\_{LA\\_Act}(t)}{\frac{\partial P\_{RA\\_Act}(t)}{\partial t} + P\_{RA\\_Act}(t)} &= F\_{RA\\_Act}(t) \cdot \begin{pmatrix} \mu\_{a2} \frac{\partial V\_{LA}(t)}{\partial t} + E\_{3LA} \\ \mu\_{a2} \frac{\partial V\_{RA}(t)}{\partial t} + E\_{3LA}(t) \end{pmatrix} \\ \tau \frac{\frac{\partial P\_{RV\\_Act}(t)}{\partial t}}{\frac{\partial T}{\partial t}} + P\_{RV\\_Act}(t) &= F\_{RV\\_Act}(t) \cdot \begin{pmatrix} \mu\_{v} \frac{\partial V\_{RV}(t)}{\partial t} + E\_{3RV} \end{pmatrix}. \end{aligned} \tag{A13}$$

Here, *P\*\*\_Act* are active parts of chamber pressures, and *E*3, μ*v*, and μ*a2* are constant parameters. *F\*\*\_Act* depended on activation time-functions *e*(t) and on the volumes *V* providing the Starling's law for the atria and the right ventricle; *kV*, *VMin*, and *VMax* are the parameters for the pressure–volume relation.

$$\begin{array}{l} F\_{LA\\_Act} = \mathfrak{e}\_{\mathfrak{a}}(t) \cdot \left( k\_{V} + (1 - k\_{V}) \cdot \left( \frac{V\_{LA} - V\_{LA\\_Minir}}{V\_{L\_{\text{M}}\text{Max}} - V\_{LA\\_Minir}} \right) \right), \\ F\_{RA\\_Act} = \mathfrak{e}\_{\mathfrak{a}}(t) \cdot \left( k\nu + (1 - k\nu) \cdot \left( \frac{V\_{RA} - V\_{RA\\_Minir}}{V\_{R\underline{a}}\text{Max} - V\_{RA\\_Minir}} \right) \right), \\ F\_{RV\\_Act} = \mathfrak{e}\_{\mathbb{P}}(t) \cdot \left( k\nu + (1 - k\nu) \cdot \left( \frac{V\_{RV\\_Mark} - V\_{RV\\_Minir}}{V\_{RV\\_Mark} - V\_{RV\\_Minir}} \right) \right). \end{array} \tag{A14}$$

Activation functions were set as follows

$$\begin{aligned} c\_d &= \left\{ \begin{array}{l} 0.5(1 - \cos\left(2\frac{(t - t\_a)\pi}{T\_a}\right) \end{array} \right\}, t\_a \le t < t\_a + T\_a\\ &0, \text{ otherwise},\\ c\_v &= \left\{ \begin{array}{l} 0.5\left(1 - \cos\left(2\frac{(t - t\_v)\pi}{T\_v}\right) \right) \cdot t\_v \le t < t\_v + T\_v\\ 0, \text{ otherwise}, \end{array} \right. \end{aligned} \tag{A15}$$

where *ta* and *tv* are the times of contraction initiations for the atria and ventricles, respectively; and *Ta, Tv* are the systole durations.

The full system of ordinary differential equations for other hemodynamic variables is

*dVRV dt* = *QiRV* − *QoRV*, *dVRA dt* = *QiRA* − *QiRV*, *dVLV dt* = *QiLV* − *QoLV* + *CLV dPLV dt* , *dVLA dt* = *QiLA* − *QiLV*, *LiLV dQiLV dt* + *RiLVQiLV* = *PLA* − *PLV*, *LoLV dQoLV dt* + *RoLVQoLV* = *PLV* − *PA*1, *LA dQA dt* + *RAQA* = *PA*1 − *PA*2, *QiRA* = *PV*−*PRA RiRA* , *QiRV* = *PRA*−*PRV RiRV* , *LoRV dQoRV dt* + *RoRVQoRV* = *PRV* − *PAPulm*,*CA*1 *dPA*1 *dt* = *QoLV* − *QA*, *CA*2 *dPA*2 *dt* = *QA* − *PA*2−*PV Rper* − *PA*2−*PC RC* , *CV dPV dt* = *PA*2−*PV Rper* − *QiRA*, *CC dPC dt* = *PA*2−*PC RC* , *CAPulm dPAPulm dt* = *QoRV* − *PAPulm*−*PVPulm RperPulm* , *CVPulm dPVPulm dt* = *QiLA* + *PAPulm*−*PVPulm RperPulm* . (A16)

*QiRV*, *QoRV*, *QiLV*, *QoLV*, *QiRA*, and *QiLA* are blood flows through the tricuspid valve, pulmonary valve, mitral valve, aortic valve, and the flows through systemic and pulmonary veins, respectively; *QA* is the arterial flow. *R* and *L* are hydraulic and inertial resistances of the vessels and chamber entrances (subindex *i*) and exits (subindex *o*); *Rper* and *RperPulm* are the peripheral vascular resistances of systemic and pulmonary circulation systems; *C* represents the compliances of the vascular reservoirs presented in the model; *RC* and *CC* characterize the viscoelastic properties of the systemic arteries. Indexes *A*1, *A*2, *V*, *APulm*, and *VPulm* corresponded to the aorta, large arteries, systemic veins, pulmonary arteries, and veins, respectively.

#### *A.4. Local Hydraulic Valve Resistances in the Cases of the Valve Pathologies*

The equations for the flows through the aortic and mitral valves in Equation (A16) were modified introducing the quadratic hydraulic valvular resistance and the dependencies of the resistances on orifice area.

$$\begin{aligned} \mathcal{L}\_{\*\*} & \stackrel{S\_{\rm S}}{\mathop{\rm S}} \frac{dQ\_{\*\*}}{dt} + \left( R\_{\*\*} + R\_{\rm s\bar{q}} |Q\_{\*\star}| \right) Q\_{\*\*} = \Delta p, \\ R\_{\rm sq} &= \mathcal{C}\_{\rm sq} \varepsilon\_{\rm Re} \mathcal{L}\_{\rm sq} S^{-2}, \\ \varepsilon\_{\rm Re} &= \sum\_{i=0}^{5} b\_{i} \left( l\_{\rm S}(\rm Re) \right)^{i}, \\ \mathcal{L}\_{\rm sq} &= \left( \mathcal{C}\_{\rm \bar{\zeta}} \Big( 1 - \frac{S}{S\_{0}} \Big)^{0.375} + 1 - \frac{S}{S\_{0}} \right)^{2}. \end{aligned} \tag{A17}$$

*S* is the orifice area, and *S*0 is an area of a completely open valve in healthy conditions; Re is the Reynolds number, *C*ζ and *bi* are parameters.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Dynamics of Periodic Waves in a Neural Field Model**

#### **Nikolai Bessonov 1, Anne Beuter 2,3, Sergei Trofimchuk 4 and Vitaly Volpert 5,6,7,\***

Received: 12 March 2020; Accepted: 29 June 2020; Published: 2 July 2020

**Abstract:** Periodic traveling waves are observed in various brain activities, including visual, motor, language, sleep, and so on. There are several neural field models describing periodic waves assuming nonlocal interaction, and possibly, inhibition, time delay or some other properties. In this work we study the influences of asymmetric connectivity functions and of time delay for symmetric connectivity functions on the emergence of periodic waves and their properties. Nonlinear wave dynamics are studied, including modulated and aperiodic waves. Multiplicity of waves for the same values of parameters is observed. External stimulation in order to restore wave propagation in a damaged tissue is discussed.

**Keywords:** neural field model; integro-differential equation; waves; brain stimulation

#### **1. Introduction**

#### *1.1. Brain Activity and Periodic Travelling Waves*

The brain displays a variety of highly nonlinear, complex dynamics across multiple spatial and temporal scales [1]. About 86 × 109 neurons of the human brain entertain complex and fluctuating interactions. Understanding the dynamics of these interactions and the mechanisms underlying their control remains a technical and theoretical challenge. In other words, when healthy brain processes evolve toward abnormal and pathological states as a result of disease, degeneration or traumatic injury, how can a therapeutic intervention be used to reposition the control parameters and guide the dynamics back toward a healthy state?  While the connectome is focused on anatomical connections, the dynamics of the networks are represented by functional connections (also called the dynome [1]). For example, brain functional connections described as a graph explore how signals are transmitted along neuroanatomical pathways and interact with local dynamics. Functional connections are often investigated via modeling [1]. One possibility is to use mathematical models to identify how an outside intervention such as neural electrical stimulation can modify local dynamics and how local dynamics will in turn affect other brain regions.

Cortical brain dynamics are investigated by means of periodic traveling waves (TW) characterized by their speed and frequency [3]. They describe the distribution of electric potential in the brain cortex. They are measured as a mean field potential (averaged macroscopic level). Propagating waves are observed during various types of brain activity. They provide subthreshold depolarization to

individual neurons and increase their spiking probability. According to Muller et al. [3], TW "travel over spatial scales that range from the mesoscopic (single cortical areas and millimetres of cortex) to the macroscopic (global patterns of activity over several centimeters) and extend over temporal scales from tens to hundreds of milliseconds." It has been proposed that TW mediate information transfer in the cortex.

Propagating waves increase the probability of neuron firing due to depolarization of neuronal membrane [4]. In [3] it was suggested that TW can be "spontaneously generated by recurrent circuits or evoked by external stimuli and travel along brain networks at multiple scales, transiently modulating spiking and excitability as they pass." The phase relations between oscillations in different cortical regions produce the TW, and depending on the distance, axonal conduction delays can reach up to tens of milliseconds. These TW correlate with the subject's performance, propagate in specific directions and synchronize distributed cortical networks that are communicating [5].

Botella-Soler et al. [6] identified for each subject the set of intracranial contacts that showed a larger percentage of detected events during slow wave cortical activity. They called these contacts "hubs" because the slow wave events in their travel through the cortical networks seemed to have a great probability of passing through the region close to the contact. Using probabilities, they were able to reconstruct a preferential propagation network for each subject. Slow waves have been reported to propagate across cortical areas at about 1m/s with multiple propagation paths and several points of origin. It seems that the slow waves have a preference to start in the prefrontal cortex and to end in posterior and temporal regions of the cortex. These waves appear to shape and strengthen neuronal networks.

Let us also note that time delays are intrinsic to the dynamics of brain networks, nodes and edges. Propagation speed along axons depends on axonal length and diameter. Conduction times along neural circuits also depend on degree of fiber myelination. The introduction of time delays in models can lead to significant changes in brain dynamics. They can be averaged or treated as distributed delays. In this paper we address the question of the effects of these delays in the dynamics of periodic cortical waves.

Thus, according to the biological observations, TW propagate in the cortex, activating and coordinating different parts of the brain. In this work, we study some of their properties. In the next section, we introduce the model. Then we present stability analysis which determines the conditions of wave appearance. Their nonlinear dynamics will be discussed in Section 3.

#### *1.2. Neural Field Model*

Neural field models were first introduced in [7]. Periodic traveling waves are described by several models (see [8–12] and Appendix A). In this work we consider one equation model with delay, and we discuss two mechanisms of the emergence of such waves, which were not sufficiently investigated previously. The first mechanism is related to the asymmetric connectivity functions [13], and the second one is determined by the delay in the response function. It is known that the loss of stability of the homogeneous in the space solution in this model does not lead to the bifurcation of periodic waves [8]. We show that they still appear for some larger values of time delay. We consider the one-dimensional neural field equation for the electric potential in the brain cortex written in the form

$$\frac{\partial \mu}{\partial t} = D \frac{\partial^2 \mu}{\partial \mathbf{x}^2} + \mathcal{W}\_d - \mathcal{W}\_i - \sigma \mathbf{u}\_\prime \tag{1}$$

where *D* is the diffusion coefficient; *Wa* and *Wi* are given by the expressions

$$\mathcal{W}\_{\mathbf{a}}(\mathbf{x},t) = \int\_{-\infty}^{\infty} \phi\_{\mathbf{a}}(\mathbf{x}-\mathbf{y}) \mathcal{S}\_{\mathbf{a}}\left(\mathbf{u}\left(\mathbf{y},t-\frac{|\mathbf{x}-\mathbf{y}|}{q\_{\mathbf{a}}}-\tau\_{\mathbf{a}}\right)\right) d\mathbf{y},\tag{2}$$

$$\mathcal{W}\_{i}(\mathbf{x},t) = \int\_{-\infty}^{\infty} \phi\_{i}(\mathbf{x}-\mathbf{y})\mathcal{S}\_{i}\left(\mu\left(y\_{i}t - \frac{|\mathbf{x}-\mathbf{y}|}{q\_{i}} - \tau\_{i}\right)\right) dy\tag{3}$$

and characterize neuron activation and inhibition. These expressions describe the intensity of signal coming from all points *y* to the point *x*; *Sa*(*u*) and *Si*(*u*) are smooth functions; *qa* and *qi* are the excitation speeds; |*x* − *y*|/*qa*,*i* is the time delay due to the excitation propagation from the point *y* to the point *x*; *φa* and *φi* are the connectivity functions,

$$\phi\_{\mathfrak{d}}(r) = \begin{cases} a\_1 e^{-b\_1 r} & , \ r > 0 \\\ a\_3 e^{b\_3 r} & , \ r < 0 \end{cases}, \quad \phi\_i(r) = \begin{cases} a\_2 e^{-b\_2 r} & , \ r > 0 \\\ a\_4 e^{b\_4 r} & , \ r < 0 \end{cases} \tag{4}$$

where *ai*, *bi* are some positive constants. Response functions *Sa*(*u*) and *Si*(*u*) are non-negative, non-decreasing functions usually considered as sigmoid-type functions; *τa* and *τi* are time delays in neuron response to the activating and inhibitory signals. The last term on the right-hand side of Equation (1) describes signal decay with a decay rate *σ* > 0.

Neural field models are often considered without the diffusion term in the studies of both brain oscillations and traveling waves (see, e.g., [14–17]). In the case of small diffusion coefficients, its influence is not essential [18]. It can describe the ephaptic effect [19], ion diffusion and gap junction. The integral terms in Equation (1) characterize nonlocal neuron communication due to axons and dendrites. The speeds *qa* and *qi* of electric impulse propagation along axons are of the order 2–4 m/s, while the speed of wave propagation is one-two orders of magnitude less [20] (p. 213). Thus, we will consider a large speed limit *qa* = *qi* → ∞:

$$\frac{\partial \mu}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + \int\_{-\infty}^{\infty} \left( \phi\_t(\mathbf{x} - \mathbf{y}) \mathbb{S}\_{\mathbf{z}} \left( u \left( y, t - \tau\_{\mathbf{z}} \right) \right) - \phi\_i(\mathbf{x} - \mathbf{y}) \mathbb{S}\_i \left( u \left( y, t - \tau\_i \right) \right) \right) dy - \sigma u. \tag{5}$$

Up to the diffusion term, which is not very essential for small diffusion coefficients, this model is a particular case of the two equation model considered in [8]. If the kernels in the two equations are the same, then the system can be reduced to the single equation.

Different neural field models describe the propagation of periodic traveling waves (see Appendix A). In this work we will consider two other mechanisms of their emergence. One of them is related to asymmetric kernels *φa* and *φi*, whose existence is confirmed by the experimental observations and used in theoretical considerations [13,21]. Another one is determined by the secondary bifurcation due to time delay for symmetric connectivity functions. Similarly to [8], we observe that the loss of stability of the homogeneous in space stationary solution leads to the appearance of periodic time oscillations independent of the space variable or of stationary periodic in space solutions. Periodic traveling waves bifurcate in the instability region, and they are unstable close to the bifurcation point. We will see that they can become stable under further change of parameters.

We studied the dynamics of the periodic waves, including their non-uniqueness for the same values of parameters. This property seems to us important, since the co-existence of waves with different frequencies is experimentally observed.

#### **2. Stability**

#### *2.1. Linearization and Eigenvalues*

In this section we will consider the Equation (5) on the interval 0 < *x* < *L* with periodic boundary conditions. We extend the function *u*(*x*, *t*) by periodicity on the whole axis, −∞ < *x* < ∞, so that the integrals in Equation (5) are well defined. Let *u*0 be a solution of the equation

$$
\phi\_a^\* S\_a(\mu) + \phi\_i^\* S\_i(\mu) - \sigma \mu = 0,
$$

where *φ*∗ *a* = ∞ −∞ *φa*(*x*)*dx*, *φ*∗ *i* = ∞ −∞ *φi*(*x*)*dx*. Then *u*0 is a stationary solution of Equation (5). Linearizing this equation about *u*0, we obtain the eigenvalue problem:

$$D\upsilon'' + S\_a'(u\_0)e^{-\lambda\tau\_0} \int\_{-\infty}^{\infty} \phi\_a(\mathbf{x} - y)\upsilon(y) dy - S\_i'(u\_0)e^{-\lambda\tau\_i} \int\_{-\infty}^{\infty} \phi\_i(\mathbf{x} - y)\upsilon(y) dy - \sigma\upsilon = \lambda\upsilon. \tag{6}$$

Here, *v* is a small perturbation of the stationary solution *u*0. Applying the Fourier transform, we get

$$S\_a'(\mathfrak{u}\_0)\varepsilon^{-\lambda\tau\_0}\tilde{\phi}\_a(\vec{\xi}) - S\_i'(\mathfrak{u}\_0)\varepsilon^{-\lambda\tau\_i}\tilde{\phi}\_i(\vec{\xi}) - D\xi^2 - \sigma = \lambda,\tag{7}$$

where

$$\begin{aligned} \tilde{\phi}\_{4}(\xi) &= \frac{a\_{1}b\_{1}}{b\_{1}^{2} + \xi^{2}} + \frac{a\_{3}b\_{3}}{b\_{3}^{2} + \xi^{2}} + i\xi \left( \frac{a\_{1}}{b\_{1}^{2} + \xi^{2}} - \frac{a\_{3}}{b\_{3}^{2} + \xi^{2}} \right), \\ \tilde{\phi}\_{i}(\xi) &= \frac{a\_{2}b\_{2}}{b\_{2}^{2} + \xi^{2}} + \frac{a\_{4}b\_{4}}{b\_{4}^{2} + \xi^{2}} + i\xi \left( \frac{a\_{2}}{b\_{2}^{2} + \xi^{2}} - \frac{a\_{4}}{b\_{4}^{2} + \xi^{2}} \right) \end{aligned}$$

are Fourier transforms of the functions *φa* and *φi*, respectively. Let us note that Fourier transform of non-integrable functions is considered in the sense of generalized functions. Using Fourier series instead of Fourier transform allows one to use a classical function instead of generalized functions. However, the advantage of the Fourier transform is that the periodicity of solutions is not imposed.

Set *λ* = *iν*. Separating the real and imaginary parts in Equation (7), we obtain:

*S a*(*u*0) cos(*ντa*) - *a*1*b*1 *b*2 1 + *ξ*2 + *a*3*b*3 *b*2 3 + *ξ*2 + *S a*(*u*0)*ξ* sin(*ντa*) - *a*1 *b*2 1 + *ξ*2 − *a*3 *b*2 3 + *ξ*2 − (8) *S i*(*u*0) cos(*ντi*) - *a*2*b*2 *b*2 2 + *ξ*2 + *a*4*b*4 *b*2 4 + *ξ*2 − *S i*(*u*0)*ξ* sin(*ντi*) - *a*2 *b*2 2 + *ξ*2 − *a*4 *b*2 4 + *ξ*2 − *Dξ*2 = *σ*, − *S a*(*u*0) sin(*ντa*) - *a*1*b*1 *b*2 1 + *ξ*2 + *a*3*b*3 *b*2 3 + *ξ*2 + *S a*(*u*0)*ξ* cos(*ντa*) - *a*1 *b*2 1 + *ξ*2 − *a*3 *b*2 3 + *ξ*2 + (9) *S i*(*u*0) sin(*ντi*) - *a*2*b*2 *b*2 2 + *ξ*2 + *a*4*b*4 *b*2 4 + *ξ*2 − *S i*(*u*0)*ξ* cos(*ντi*) - *a*2 *b*2 2 + *ξ*2 − *a*4 *b*2 4 + *ξ*2 = *ν*.

#### *2.2. Symmetric Connectivity Functions with Time Delay*

If the symmetry condition

$$a\_1 = a\_3, a\_2 = a\_4, b\_1 = b\_3, b\_2 = b\_4 \tag{10}$$

is satisfied, then Equations (8) and (9) are as follows:

$$S\_a'(u\_0)\cos(\nu\tau\_d)\frac{a\_1b\_1}{b\_1^2+\zeta^2} - S\_i'(u\_0)\cos(\nu\tau\_i)\frac{a\_2b\_2}{b\_2^2+\zeta^2} - D\zeta^2/2 = \sigma/2,\tag{11}$$

$$-S\_a'(\mu\_0)\sin(\nu\tau\_d)\frac{a\_1b\_1}{b\_1^2+\xi^2} + S\_i'(\mu\_0)\sin(\nu\tau\_i)\frac{a\_2b\_2}{b\_2^2+\xi^2} = \nu/2. \tag{12}$$

We can express *ξ*2 through *ν* from the last equation and substitute them into Equation (11). The resulting equation with respect to *ν* should be solved numerically or asymptotically. Since the calculations are sufficiently complex, we will consider here a simplified case where *τa* = 0. Numerical simulations show that periodic waves can exist in this case (Section 3). Our aim here is to analyze their bifurcations from the homogeneous in space solution. Assuming that *ν* ≥ 0, we also get the conjugate solution −*ν*.

If *τa* = 0, then Equations (11) and (12) are as follows:

$$S\_a'(u\_0)\frac{a\_1 b\_1}{b\_1^2 + \zeta^2} - S\_i'(u\_0)\cos(\nu \tau\_i) \frac{a\_2 b\_2}{b\_2^2 + \zeta^2} - D\zeta^2/2 = \sigma/2,\tag{13}$$

$$\frac{1}{2}S\_i'(\mu\_0)\sin(\nu\tau\_i)\frac{a\_2b\_2}{b\_2^2+\xi^2} = \nu/2. \tag{14}$$

By the stability boundary we will understand here the value *σ* = *σ*0 such that system (13), (14) does not have solutions for *σ* > *σ*0, and it has at least one solution for *σ* < *σ*0. This stability boundary *σ*0 is uniquely defined and it can be negative.

#### **Proposition 1.** *At the stability boundary, ξ* = 0 *or ν* = 0*.*

**Proof.** Consider Equation (14) as an equation with respect to *ν* for a fixed *ξ*. It has a solution *ν* = 0 if and only if

$$2\pi\_i S\_i'(u\_0) \frac{a\_2 b\_2}{b\_2^2 + \xi^2} > 1. \tag{15}$$

Suppose that this condition is satisfied for some *ξ*0 > 0, and denote by *ν*(*ξ*) solution of Equation (14) in the vicinity of *ξ* = *ξ*0 > 0, and *ν*0 = *ν*(*ξ*0) > 0. Assume that the corresponding value *σ*0 belongs to the stability boundary.

Taking into account Equation (14), we can write Equation (13) as follows:

$$2S\_a'(\mu\_0)\frac{a\_1b\_1}{b\_1^2+\xi^2} - \nu \cot(\nu \tau\_i) - D\xi^2 = \sigma. \tag{16}$$

Let *ξ*1 be sufficiently close to *ξ*0 such that condition (15) remains satisfied, and *ξ*1 < *ξ*0. Considering *ν*(*ξ*0) < *π*/*τi*, we get *ν*(*ξ*1) > *ν*(*ξ*0). Therefore,

$$\nu(\boldsymbol{\xi}\_1)\cot(\nu(\boldsymbol{\xi}\_1)\boldsymbol{\tau}\_i) < \nu(\boldsymbol{\xi}\_0)\cot(\nu(\boldsymbol{\xi}\_0)\boldsymbol{\tau}\_i) \dots$$

Set

$$F(\xi) = 2S\_a'(\mu\_0) \frac{a\_1 b\_1}{b\_1^2 + \xi^2} - \nu(\xi) \cot(\nu(\xi) \tau\_i) - D\xi^2.$$

Thus, *F*(*ξ*1) > *F*(*ξ*0). Hence, *σ*0 = *F*(*ξ*0) cannot belong to the stability boundary since *σ*1 = *F*(*ξ*1) > *σ*0, and system (13), (14) has a solution for *σ* = *σ*1. This contradiction shows that *ξ*0 = 0 when *ν*0 = 0.

Let 2*π*/*τi* < *ν*(*ξ*0) < 5*π*/(2*τi*). Then there is another solution *ν*1(*ξ*0) of Equation (14) such that *π*/(2*τi*) < *ν*1(*ξ*0) < *π*/*τi*. Therefore, cos(*ν*1(*ξ*0)*τi*) < 0, cos(*ν*(*ξ*0)*τi*) > 0, and

$$
\omega\_0 = F(\xi\_0) < 2S\_a'(\mu\_0) \frac{a\_1 b\_1}{b\_1^2 + \xi\_0^2} - \nu\_1(\xi\_0) \cot(\nu\_1(\xi\_0)\tau\_i) - D\xi\_0^2.
$$

Hence, *σ*0 cannot belong to the stability boundary.

Finally, in all other cases with *ν*(*ξ*0) > 5*π*/(2*τi*) we obtain a contradiction similarly to the two cases considered before. This contradiction shows that at the stability boundary, either *ξ* = 0 or *ν* = 0.

Let us now determine the stability boundary with respect to time delay for a fixed *σ*. If *ν* = 0, then from Equation (13) we get

$$
\sigma = \frac{a\_1}{b\_1^2 + \xi^2} - \frac{a\_2}{b\_2^2 + \xi^2} - D\xi^2,\tag{17}
$$

where *α*1 = 2*S a*(*u*0)*a*1*b*1, *α*2 = 2*S i* (*u*0)*a*2*b*2. Assuming that *ν* = 0, we set *βi* = *αi*/(*b*2 *i* + *ξ*2), *i* = 1, 2. Then *sin*(*ντi*) = *ν*/*β*2,

$$
\beta\_1 - \sqrt{\beta\_2^2 - \nu^2} = \sigma, \quad \nu^2 = \beta\_2^2 - (\beta\_1 - \sigma)^2.
$$

**Example 1.** *Consider the values of parameters: S a*(*u*0) = *S i* (*u*0) = 20, *a*1 = *a*2 = 4, *b*1 = 40, *b*2 = 20, *α*1 = 6400, *α*2 = 3200*, σ* = 0.01*, Then we find*

$$\beta\_1 = 0, \beta\_1 = 4, \beta\_2 = 8, \nu^2 \approx 48, \nu \approx 6.93, \ \tau\_i = \frac{1}{\nu} \arcsin \frac{\nu}{\beta\_2} = 0.151,$$

*and*

$$
\mu\_1^\pi = \pi, \beta\_1 = 3.975, \beta\_2 = 7.805, \nu^2 \approx 45.123, \nu \approx 6.717, \quad \tau\_i = \frac{1}{\nu} \arcsin \frac{\nu}{\beta\_2} = 0.154
$$

*Hence, the homogeneous in space oscillations appear for a smaller value of time delay than periodic traveling waves. Therefore, these waves are unstable in the vicinity of the bifurcation point. We will see in the next section that they become stable for larger values of time delay.*

#### *2.3. Asymmetric Connectivity Function without Time Delay*

If conditions (10) do not not hold, then eigenvalue (7) has a nonzero imaginary part, *λ*(*ξ*) = *α*(*ξ*) + *iν*(*ξ*),

$$\alpha(\xi) = \frac{s\_a a\_1 b\_1}{b\_1^2 + \xi^2} + \frac{s\_a a\_3 b\_3}{b\_3^2 + \xi^2} - \frac{s\_i a\_2 b\_2}{b\_2^2 + \xi^2} - \frac{s\_i a\_4 b\_4}{b\_4^2 + \xi^2} - D\xi^2 - \sigma,$$

$$\nu(\xi) = s\_a \xi \left(\frac{a\_1}{b\_1^2 + \xi^2} - \frac{a\_3}{b\_3^2 + \xi^2}\right) - s\_i \xi \left(\frac{a\_2}{b\_2^2 + \xi^2} - \frac{a\_4}{b\_4^2 + \xi^2}\right),\tag{18}$$

*sa* = *S a*(*u*0),*si* = *S i* (*u*0). Condition *α*(*ξ*) = 0 determines the stability boundary,

$$\sigma = \frac{s\_4 a\_1 b\_1}{b\_1^2 + \xi^2} + \frac{s\_4 a\_3 b\_3}{b\_3^2 + \xi^2} - \frac{s\_i a\_2 b\_2}{b\_2^2 + \xi^2} - \frac{s\_i a\_4 b\_4}{b\_4^2 + \xi^2} - D\xi^2 \equiv \Phi\_2(\xi). \tag{19}$$

If *σ* < Φ2(*ξ*) for some values of *ξ*, then the solution loses its stability due to a pair of complex conjugate eigenvalues *λ*(*ξ*) = *α*(*ξ*) ± *iν*(*ξ*). For *α* = 0 corresponding to the stability boundary, the bounded solution of the linearized equation can be written as follows:

$$u(\mathfrak{x},t) = e^{i\nu t}e^{i\mathfrak{j}\mathfrak{x}} + e^{-i\nu t}e^{-i\mathfrak{j}\mathfrak{x}} = \cos(\nu t + \mathfrak{j}\mathfrak{x}).$$

We find *ξ* from equality *α*(*ξ*) = 0 and *ν* from (18). The frequency of the periodic wave equals *ξ* and its speed *c* = −*ν*/*ξ* can be determined from (18):

$$c = -s\_d \left(\frac{a\_1}{b\_1^2 + \xi^2} - \frac{a\_3}{b\_3^2 + \xi^2}\right) + s\_l \left(\frac{a\_2}{b\_2^2 + \xi^2} - \frac{a\_4}{b\_4^2 + \xi^2}\right). \tag{20}$$

Different waves with the frequencies satisfying condition *σ* < Φ(*ξ*) can exist for the same values of parameters. Their speed can be increasing or decreasing functions of frequency according to (20).

#### **3. Numerical Results**

Numerical simulations of Equation (5) in a bounded interval 0 ≤ *x* ≤ *L* with periodic boundary conditions will be started in the case without time delay and continued with the case of time delay in the response functions *Sa*(*u*) = *Sb*(*u*) = arctan(*hu*), *h* > 0. We will finish this section with modeling of the stimulation of the damaged tissue in order to restore wave propagation.

The numerical method uses an implicit finite difference scheme with Thomas's algorithm for the inversion of the tridiagonal matrix. Initial conditions are considered either in the form of a piece-wise constant function or a sinus-like function with a given periodicity. The continuation method is used to follow the branches of solutions. In this case, the value of some parameter gradually changes. As an initial condition for the new value of parameter, we take the result of the simulation obtained for the previous value of parameter. The time step was usually taken as 0.05 and the space step 0.005. The accuracy of numerical simulations was controlled by decreasing the time and space steps. The typical values of parameters used in numerical simulations are presented in Table A1. Their specific values are give below.

#### *3.1. Wave Propagation without Time Delay*

In the case without time delay (*τa* = *τi* = 0) and with symmetric connectivity functions *φa* and *φi*, where *a*1 = *a*3, *a*2 = *a*4, *b*1 = *b*3, *b*2 = *b*4, periodic in space stationary solutions bifurcate from the constant solution. Linear stability analysis shows that traveling waves with nonzero speed do not bifurcate in this case.

If the connectivity functions are not symmetric, then traveling waves with nonzero speed are observed (see Figures 1 and 2). Linear stability analysis allows the determination of the wave speed and frequency near the bifurcation point. Let us consider an example with the following values of parameters: *a*1 = 3, *a*2 = 4, *a*3 = 3, *a*4 = 1, *b*1 = 40, *b*2 = 20, *b*3 = 40, *b*4 = 20, *D* = 10−3, *h* = 20, *L* = 2. Analysis of the function *α*(*ξ*) (Section 2.3) shows that it has a maximum at *ξ* = 20.53. It is positive for *σ* < *σc* = 0.4 resulting in the bifurcation of the periodic wave with the speed given by Equation (20). For the chosen values of parameters, we get *c* = −0.54. Numerical simulations show the emergence of a periodic wave with 7 periods in the interval [0, *L*] and with the corresponding frequency *ξ* = 21.98. Since the number of periods is an integer, the frequency of the wave observed in numerical simulations is not precisely equal to the analytical value. The wave with the closest frequency emerges since the corresponding eigenvalue has the maximal real part in comparison with the waves with other periods. The speed of this wave *c* = −0.52 is close to the analytical value.

For the values of parameters far from the stability boundary, two types of solutions are observed in numerical simulations: periodic waves with a constant speed (Figure 1) and aperiodic waves with oscillating speed (Figure 2). In the first case, the wave has conventional form *w*(*x* − *ct*), where *w*(*x*) is a periodic in space function and *c* is a constant. In the second case, the amplitude of spatial peaks and the wave speed oscillate.

**Figure 1.** A snapshot of periodic wave described by Equation (5) for *τa* = *τi* = 0 (left) and a zoom-in on the lower part of the graph (right). The values of parameters are as follows: *D* = 10−4, *σ* = 0.01, *L* = 2, *a*1 = *a*2 = 0.6, *a*3 = *a*4 = 4, *b*1 = *b*3 = 40, *b*2 = *b*4 = 20. Here *Sa*(*u*) = *Si*(*u*) = arctan(*hu*), and *h* = 20. Here and in all examples below, the stationary solution *u*0 = 0.

**Figure 2.** A snapshot of aperiodic wave described by Equation (5) for *τa* = *τi* = 0 (left) and a zoom-in on the lower part of the graph (right). Periodic (Figure 1) and aperiodic waves can exist for the same values of parameters. Convergence of solutions to one of them is determined by the initial conditions. The values of parameters are as follows: *D* = 10−4, *σ* = 0.01, *L* = 2, *a*1 = *a*2 = 0.6, *a*3 = *a*4 = 4, *b*1 = *b*3 = 40, *b*2 = *b*4 = 20. Here *Sa*(*u*) = *Si*(*u*) = arctan(*hu*), and *h* = 20.

Different periodic regimes can co-exist for the same values of parameters. Figure 3 shows the dependence of the wave speed on the values *a*1(= *a*2) for all other parameters fixed. There are three branches of solutions corresponding to different spatial frequencies. Thus, there are three different periodic waves for the same values of parameters with the speeds depending on their frequency.

**Figure 3.** The speed of periodic waves for different values of the parameters *a*1 = *a*2. Three curves correspond to different values of the spatial frequency: (1) solution with 12 periods in the interval [0, *L*], (2) 13 periods, (3) 14 periods. Branching in curve 1 shows the transition to modulated oscillations with the maximal and minimal values of the oscillating speed. The values of parameters are as follows: *D* = 10−4, *σ* = 0.01, *L* = 2, *a*3 = *a*4 = 4, *b*1 = *b*3 = 40, *b*2 = *b*4 = 20. Here *Sa*(*u*) = *Si*(*u*) = arctan(*hu*), and *h* = 20.

For a sufficiently small *a*1, transition to modulated waves can occur (curve 1) with the amplitude and speed depending on time. Such solutions can be qualitatively approximated by the function *u*(*x*, *t*)=(*k*1 + sin(*k*2*x* + *k*3*t*)) sin(*k*4*x* + *k*5*t*), where *ki* are some constants. For even smaller values of

*a*1, transition to aperiodic oscillations is observed for all three branches of solutions. These oscillations can coexist with periodic or modulated periodic waves for the same values of parameters.

#### *3.2. Time Delay and Symmetric Connectivity Functions*

#### 3.2.1. Initial Conditions

As it was discussed in the previous section, the loss of stability of the homogeneous in space stationary solution leads either to the homogeneous in space time oscillations or to the stationary periodic in space solutions. Consider time delay as a bifurcation parameter. If it exceeds a critical value, then time oscillations emerge. Periodic traveling waves bifurcate for a larger value of time delay, and they are unstable in the vicinity of the bifurcation point. Numerical simulations show that they can become stable under a further increase of *τ*. Nevertheless, periodic time oscillations homogeneous in space remain stable. Therefore, we need to choose some particular initial conditions in order to get periodic traveling waves.

The simulations presented in this section were carried out with two types of initial conditions. In the first case, we consider the equation

$$\frac{\partial \mu}{\partial t} = D \frac{\partial^2 \mu}{\partial x^2} + I(\mathbf{x}, t) \tag{21}$$

on the time interval 0 ≤ *t* ≤ *T*0. Here *I*(*x*, *t*) = *I*0 cos(*px* + *qt*), and the initial condition *u*(*x*, 0) = 0. The result of this simulation is considered as the initial condition for Equation (5). We set *I*0 = 0.5, *q* = 0.015, *T*0 = 20; the value of *p* is taken 3, 6, 9 depending on the required space periodicity.

The second type of initial conditions is used in the continuation method. The results of the simulations of Equation (5) for some values of parameters are used as initial conditions for some other values of parameters.

#### 3.2.2. Multiplicity of Waves and Parameter Dependence

An example of periodic traveling waves with different initial conditions and the same values of parameters are shown in Figure 4. We observe one, two and three-period waves generated by function *I*(*x*, *t*) for *p* = 3, 6 and 9. These waves have different speeds and amplitudes. The waves with larger wavelengths have higher speeds and amplitudes. If time delay is sufficiently small, then the waves become unstable, and the transition to periodic time oscillations independent of the space variable is observed. If the value *I*0 is small enough, then the solution converges to the homogeneous time oscillations.

**Figure 4.** Three types of periodic waves observed for the same values of parameters and having different spatial frequency and speed: one-period wave with speed −0.027 (left); two-period wave with speed −0.012 (middle); three-period wave with speed −0.0094. The values of parameters are as follows: *D* = 10−4, *σ* = 0.01, *L* = 2, *τa* = 0, *τi* = 12, *a*1 = *a*2 = *a*3 = *a*4 = 4, *b*1 = *b*3 = 40, *b*2 = *b*4 = 20. Here *Sa*(*u*) = *Si*(*u*) = arctan(*hu*), and *h* = 20.

Increase of time delay leads to the increase of the wave amplitude and to the decrease of its speed (Figure 5). Wave propagation is determined by the transmission of activating signal between the neurons. Since time delay retards this transmission, the speed decreases. This result was previously obtained analytically for monotone waves [22]. Since the activating signal propagates slower with time delay, the local response and the wave amplitude increase.

**Figure 5.** The absolute value of the speed (1a, 2a) and the amplitude (1b, 2b) for the one-period wave (1a, 1b) and two-period wave (2a, 2b) depending on *τi*. Connected points correspond to stable solutions; separate points to unstable solutions. The latter lead to the appearance of stationary solutions periodic in space. The values of parameters are as follows: *D* = 10−4, *σ* = 0.01, *τa* = 0, *a*1 = *a*2 = *a*3 = *a*4 = 4, *b*1 = *b*3 = 40, *b*2 = *b*4 = 20. Here *Sa*(*u*) = *Si*(*u*) = arctan(*hu*), and *h* = 20.

Let us note that the critical value of time delay *τi* ≈ 1.5 found in numerical simulations and the value of the wave speed near the bifurcation point *c* ≈ −2 correspond to the analytical values determined in the example in Section 2.2 (Figure 5, curves 1a, 1b). The analytical value of the speed *c* = −*ν*/*ξ* = −2.13 at the bifurcation point is slightly different from the numerical value. Since the periodic wave is unstable near the bifurcation point, we can determine its speed in numerical simulations only at some distance from the critical value *τi*.

Figure 6 shows the dependence of periodic waves on the value *b*2 (equal to *b*4). The wave amplitude and speed decrease with the increase of *b*2. There is a critical value *b*2 ≈ 40 for which the speed becomes zero, and a transition to another branch of solutions is observed. These are stationary, periodic in space solutions with growing amplitude as *b*2 increases. It is interesting to note the existence of weakly oscillating time periodic solutions in a narrow interval between traveling waves and stationary solutions.

The results presented above were obtained for a single delay *τi* in the inhibition term while *τa* = 0. If we fix *τi* = 1 and increase *τa* beginning from *τa* = 0, then the amplitude and the speed of traveling waves are not monotonous. The former first decreases, passes through the minimum and then increases, while the latter increases in the beginning and decreases for larger values of the delay (not shown).

**Figure 6.** Dependence of the amplitude of the observed regimes on *b*2(= *b*4). Periodic traveling waves are located on branch 1; stationary solutions on branch 3. A short branch 2 contains weakly oscillating solutions. The values of parameters are as follows: *D* = 10−4, *σ* = 0.01, *L* = 2, *τa* = 0, *τi* = 1, *a*1 = *a*2 = *a*3 = *a*4 = 4, *b*1 = *b*3 = 40. Here *Sa*(*u*) = *Si*(*u*) = arctan(*hu*), and *h* = 20.

#### *3.3. Stimulation*

Wave propagation can be different in a damaged tissue [23]. In this section we discuss its possible restoration by external stimulation. This question was considered in [24] for monotonous wave fronts.

#### 3.3.1. Exact Solution of the Stimulation Problem

Let us write Equation (5) for the normal tissue

$$\frac{\partial \mu}{\partial t} = D \frac{\partial^2 \mu}{\partial \mathbf{x}^2} + f(\mu) - \sigma \mathbf{u}\_\prime \tag{22}$$

$$\text{where}\\
\qquad J(u) = \int\_{-\infty}^{\infty} \phi\_d(\mathbf{x} - y) \mathbf{S}\_d \left( u \left( y, t - \mathbf{r}\_d \right) \right) dy - \int\_{-\infty}^{\infty} \phi\_i(\mathbf{x} - y) \mathbf{S}\_i \left( u \left( y, t - \mathbf{r}\_i \right) \right) dy,$$

and a similar equation for the damaged tissue

$$\frac{\partial \upsilon}{\partial t} = D \frac{\partial^2 \upsilon}{\partial x^2} + f^\*(\upsilon) - \sigma \upsilon,\tag{23}$$

$$\text{where}\\
\qquad J^\*(\boldsymbol{\upsilon}) = \int\_{-\infty}^{\infty} \phi\_a^\*(\mathbf{x} - \boldsymbol{y}) S\_a^\* \left(\boldsymbol{\upsilon} \left(\boldsymbol{y}, t - \tau\_a\right)\right) d\boldsymbol{y} - \int\_{-\infty}^{\infty} \phi\_i^\*(\mathbf{x} - \boldsymbol{y}) S\_i^\* \left(\boldsymbol{\upsilon} \left(\boldsymbol{y}, t - \tau\_i\right)\right) d\boldsymbol{y}, \quad i = 1, \ldots, n$$

and ∗ denotes the functions for the damaged tissue. In the presence of stimulation *I*(*x*, *t*) this equation is as follows:

$$\frac{\partial z}{\partial t} = D \frac{\partial^2 z}{\partial x^2} + f^\*(z) - \sigma z + I(x, t). \tag{24}$$

We will choose this function such that solution *z*(*x*, *t*) of Equation (24) becomes equal solution *u*(*x*, *t*) of Equation (22). We set *z*(*x*, *t*) = *u*(*x*, *t*), and substitute *u*(*x*, *t*) in Equation (24). Then

$$I(\mathbf{x}, t) = \frac{\partial \mu}{\partial t} - D \frac{\partial^2 \mu}{\partial \mathbf{x}^2} - J^\*(\mu) + \sigma \mu = J(\mu) - J^\*(\mu).$$

Thus, the stimulation function

$$I(\mathbf{x}, \mathbf{t}) = J(\boldsymbol{\mu}) - J^\*(\boldsymbol{\mu})$$

gives solution of the complete reconstruction problem.

#### 3.3.2. Approximate Solution of the Stimulation Problem

The stimulation function suggested above uses the solution *u*(*x*, *t*) which may not be known for the patient with the damaged brain tissue. In this case, some approximate solutions of the stimulation problem can be used. Consider the integral *J*∗(*v*) in the following form:

$$J^\*(\upsilon) = \int\_{-\infty}^{\infty} \mathcal{W}(\mathbf{x}) \mathcal{W}(\mathbf{y}) \left( \phi\_d(\mathbf{x} - \mathbf{y}) \mathbb{S}\_a \left( \upsilon \left( \mathbf{y}, t - \tau\_d \right) \right) - \phi\_i(\mathbf{x} - \mathbf{y}) \mathbb{S}\_i \left( \upsilon \left( \mathbf{y}, t - \tau\_i \right) \right) \right) d\mathbf{y} \Big|$$

where *W*(*x*) = *w*0 < 1 for *x*1 ≤ *x* ≤ *x*2, *x*1, *x*2 ∈ (0, *L*), and *W*(*x*) = 1 outside the interval [*x*1, *x*2]. Hence, the connectivity function decreases if *x* or *y* belongs to the damaged area [*x*1, *x*2]. If we set *w*0 = 0, then the connectivity function vanishes if one of the two points *x* or *y* (neurons) belongs to the damaged area.

Without damage, the periodic traveling wave solution of Equation (22) and the integral *J*(*u*) can be approximated by a cosine function (Figure 7). Hence, we will look for the stimulation function in the form approximating the periodic wave: *I*(*x*, *t*) = *I*0(*x*) cos(*px* + *qt*). We set *p* = 6, *q* = 1 to approximate the frequency and the speed of the wave in the normal tissue. We set *I*0(*x*) = *i*0 for *x*1 ≤ *x* ≤ *x*2, and *I*0(*x*) = *i*1 outside the interval [*x*1, *x*2]. The results of numerical simulations are presented in Figure 7 with the comparison of the normal tissue, damaged tissue without stimulation and damaged tissue with stimulation. We observe a periodic wave in the normal tissue. The behavior of the solution is completely different in the damaged tissue. The solution is close to 0 at the damaged interval (green in the middle figure), and it oscillates periodically in time from both sides of this interval. Stimulation restores the wave propagation with the same frequency and speed. The choice of the stimulation amplitude *i*0 = 0.6 and *i*1 = 0.1 (for the example in the figure) is important. If we set *i*0 = 0.6, *i*1 = 0, the stimulation is not successful; there is no wave propagation. Thus, stimulation should be done not only inside the damaged interval but also around it.

**Figure 7.** Periodic traveling wave in the normal tissue and the integral *J*(*u*) (**left**). Snapshots of solutions in the damaged tissue for two different moments of time (**a**) and (**b**) (**middle**). Solution with stimulation becomes close to the periodic wave (**right**). The values of parameters are as follows: *D* = 10−4, *σ* = 0.01, *L* = 2, *τa* = 0, *τi* = 1, *a*1 = *a*2 = *a*3 = *a*4 = 4, *b*1 = *b*3 = 40, *b*2 = *b*4 = 20, *Sa*(*u*) = *Si*(*u*) = arctan(*hu*), *h* = 20, *p* = 6, *q* = 1, *i*0 = 0.6, *i*1 = 0.1. Green interval [0.5, 1.07] shows the damaged tissue, *w*0 = 0.

#### **4. Discussion**

Brain functioning is determined by large scale networks of epicenters (hubs) located in the cortex and connected by white matter fiber tracks. The structure of these networks depends on the particular type of the brain activity (motor, language and so on) and on the inter-individual variation. The epicenters exchange information by means of signaling along the cortex or along the white matter fibers. Apparently, this signaling occurs in the form of traveling waves. Periodic traveling waves are observed in thalamus, visual cortex, hippocampus and other parts of the brain. There is more and more evidence that they play the key roles in brain functioning. However, their exact roles and organization are not known. Their characteristics, such as speed, frequency and amplitude, can vary in wide limits (Appendix B), and it is not clear how they are initiated and stop, and how they are related to some particular types of activities. For example, periodic waves are observed in the beginning of produced speech, and they inverse their direction at the end [25]. The corresponding mechanisms governing these activities are not yet understood. These investigations are at the stage of accumulation of biological data and of elaboration of different models whose role and utility will become clearer with time.

Neural field models are widely used to study brain patterns, including stationary structures, pulses, wave fronts, and periodic waves. There are a number of models which describe periodic waves (Appendix A) since it is relatively easy to find these waves on the basis of linear stability analysis of the homogeneous in space stationary solution. The mechanisms of this instability can be related to a combination of inhibition, time delay, refractoriness (also a variant of time delay) and nonlocal interaction. In this work we study the influence of asymmetry of connectivity functions and of time delay in neuron response.

Neuron connectivity is provided by numerous axons whose density decrease approximately exponentially as a function of distance [26,27]. In modeling, the connectivity functions are usually considered to be symmetric; that is, axon connections between points *x* and *y* have the same density as between *y* and *x*. There is some evidence that connectivity can be asymmetric [13].

In this work we study periodic waves emerging due to the asymmetry of the connectivity function and due to time delay in the response function. Both of them correspond to the existing biological mechanisms. They allow us to study periodic waves in the minimal model which consists of a single integro-differential equation. Other models contain two equations (Appendix A). In the case of asymmetric connectivity function, periodic waves bifurcate due to the loss of stability of the homogeneous in space stationary solution. They become unstable, while bifurcating periodic waves are stable. This is different in the case of time delay. Periodic waves bifurcate from the unstable stationary solution because of space-independent oscillations. Increasing time delay leads to their stabilization.

Though the model contains quite many parameters, using some combination of parameters reduces their number. Next, there are some additional relations between the parameters which determine the stability boundary. Furthermore, there are some physical estimates, such as the decay rate of neuron connectivity function or the value of time delay. Finally, the wave speed and frequency, which depend on parameters, should be in some experimentally observed range. Altogether, these constraints determine some limited intervals of parameter variation (Appendix B).

#### *Cortex Damage and Stimulation*

The network of hubs related to some brain function (connectome) can be damaged because of stroke or other factors leading to partial or complete loss of the corresponding function. Post-stroke patient recovery is often incomplete, and usually limited to six months after the accident. Various stimulation techniques are discussed in the literature, but their results are controversial [28].

Post-stroke brain damage can influence propagation of brain waves between the hubs of the connectome. Traveling waves in the cortex can have several functions, including activation of some of its parts. This activation facilitates firing of individual neurons [4]. TW reflect information originating from right and left hemispheres, traveling short and long distances, and containing various time delays. In other words, TW coordinate multiple faster and slower speech events by preparing the arrival of signals traveling along white matter tracts to specific hubs. In particular, direct brain recording (ECoG) showed the presence of TW during consonant-vowel syllables' pronunciation [25]. Summarizing the results by Gross et al. [29] examining how brain waves help us make sense of speech in healthy subjects, Weaver [30] indicated the presence of different time scales in speech production from tens of milliseconds (phoneme) to hundreds of milliseconds (intonation). Today, the language networks have been identified with precision, including phonological processing, speech planning, language semantics, spatial cognition and other functions [31].

We address in this work the question about possible restoration of cortex waves by external stimulation. We show that appropriate choice of injected current allows the recovery of wave speed and frequency, possibly leading to a better communication between the hubs of the connectome. This proof of concept is the first step on the long road to a possible application of this approach to patient rehabilitation.

**Author Contributions:** Conceptualization, V.V.; Data curation, A.B.; Formal analysis, S.T.; Methodology, A.B. and V.V.; Software, N.B.; Writing original draft, V.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** A. Beuter was funded by CorStim SAS. S. Trofimchuk was partially supported by FONDECYT (Chile), project 1190712. This research was also funded by "RUDN University Program 5-100" and the French-Russian project PRC2307.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Periodic Waves in Different Models**

#### **Two equations with time delay.**

In the work [8] a single neuron population model (excitation)

$$
\pi \frac{\partial u}{\partial t} = \int\_{-\infty}^{\infty} P(\mathbf{x} - y) \Psi(\mathbf{u}(y, t - d)) d\mathbf{u} - \mathbf{u} \tag{A1}
$$

and a two population model (excitation and inhibition)

$$\tau \frac{\partial u}{\partial t} = \int\_{-\infty}^{\infty} \left( P\_{11}(\mathbf{x} - \mathbf{y}) \psi\_1(u(y, t - d)) - P\_{12}(\mathbf{x} - \mathbf{y}) \psi\_2(v(y, t - d)) \right) dy - u \tag{A2}$$

$$
\pi \frac{\partial \upsilon}{\partial t} = \int\_{-\infty}^{\infty} \left( P\_{21}(\mathbf{x} - \mathbf{y}) \psi\_1(\mathbf{u}(y, t - d)) - P\_{22}(\mathbf{x} - \mathbf{y}) \psi\_2(\upsilon(y, t - d)) \right) dy - \upsilon \tag{A3}
$$

are considered. Here *P*(*x*) and *Pij*(*x*) are symmetric positive functions. A particular example of step-wise constant functions is studied. Periodic traveling waves cannot exist for the first model. They are observed for the second model. The existence of pulse solutions in a similar model without delay was studied in [32].

#### **One equation with distributed speed and delay.**

Equation with a distributed propagation speed and time delay is considered in [9]:

$$L\left(\frac{\partial \mu}{\partial t}\right) = \mathfrak{a} \int\_0^\infty \mathfrak{g}(\upsilon) \int\_{-\infty}^\infty K(z) \mathcal{S}(\mathfrak{u}(\mathfrak{x} + z, t - |z|/\upsilon)) dz d\upsilon +$$

$$\mathfrak{F} \int\_0^\infty f(\tau) \int\_{-\infty}^\infty F(z) \mathcal{S}(\mathfrak{u}(\mathfrak{x} + z, t - \tau)) dz d\tau,\tag{A4}$$

where *L* is a second-order differential operator; the functions *K*(*z*) and *F*(*z*) include both activatory and inhibitory kernels. Different regimes are observed: periodic in time and independent of space, stationary periodic in space, periodic traveling waves.

The model with distributed time delay is considered in [10]

$$u(\mathbf{x},t) = \int\_{-\infty}^{\infty} w(\mathbf{x}-\mathbf{y})d\mathbf{y} \int\_{-\infty}^{t} \eta(t-s)f(u(\mathbf{y},\mathbf{s}-|\mathbf{x}-\mathbf{y}|/v))ds.$$

The same regimes as above and oscillating Turing structures are observed.

#### **Neural field model with linear adaptation.**

It is a two equation model

$$
\pi \frac{\partial u}{\partial t} = -u - \beta v + \int\_D w(\mathbf{x} - y) \mathbb{F}(u(y, t)) dy,\tag{A5}
$$

$$\frac{1}{a}\frac{\partial v}{\partial t} = u - v,\tag{A6}$$

where the second variable represents a linear adaptation (see [11] and the references therein). A large variety of waves and patterns are observed, including stationary periodic in space solutions, traveling waves, modulated traveling waves and stationary and oscillating bumps.

#### **Neural field model with refractoriness.**

Periodic traveling waves are also found in one-equation model without inhibition term but with neuron refractoriness (time delay after firing) [12]:

$$\frac{1}{r}\frac{\partial u}{\partial t} = -u + \left(1 - \int\_{t-1}^{t} u(x,s)ds\right)f(w\otimes u).$$

Here ⊗ denotes spatial convolution.

#### **Appendix B. The Values of Parameters**

Periodic traveling waves behave as cos(*βt* + *ξx*) with the time frequency *β*, space frequency *ξ* and speed *c* = −*β*/*ξ*. With the interval length *L* = 2 cm, we get *ξ* = 2*π*/*L* ≈ 3 cm−1.

Let us consider the example of the simulations in Figure 5 with *τ* = 1 and wave speed 0.3. If this value of *τ* corresponds to the characteristic time delay 10 ms ([9,33] Chapter 2.3.1), that is, the time unit in the simulation corresponds to 0.01 s, then *c* = 0.3 × 100 = 30 cm/s.

The time frequency *β* = *cξ* = 90 s−1 belongs to the upper limit of the observed range. The value of the wave speed, and respectively, the time frequency linearly dependent on it can be decreased by the variation of parameters *ai* and *bi*.

Connectivity functions can be estimated from the data in [26,27]. It exponentially decreases with the rate of decrease in the interval 3–10 times at the at the distance 0.03 cm. This corresponds to the exponential exp(−*μx*) with *μ* in the range 30 ÷ 40 cm−1.

Propagation speed measures vary depending of the methodology used. When macroscopic waves are recorded from EEG or from ECoG which have low spatial and high temporal resolutions, the propagation speeds varies between 1 and 10 m/s. As indicated by Muller et al. [3], these results are compatible with the range of axonal conduction speeds of myelinated white matter fibers in the cortex. However, when measuring mesoscopic waves' propagation speed using local field potential (LFP) from multielectrode arrays (MEAs) or from optical imaging signals recorded with voltage-sensitive dyes (VSDs) having high spatial and temporal resolution, the propagation speeds varies from 0.1 to 0.8 m per second, consistent with the axonal conduction speed of the unmyelinated long-range horizontal fibers within the superficial layers of the cortex, as indicated by Muller et al. [3].

#### **Summary of parameters.**

Parameters of the model are summarized in the following table. Let us note that all parameters are dimensionless. They can be recalculated to dimensional values as in the example above.

**Table A1.** Typical values of parameters of the model.

#### **References**

c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Mathematical Modeling Shows That the Response of a Solid Tumor to Antiangiogenic Therapy Depends on the Type of Growth**

#### **Maxim Kuznetsov 1,2**

Received: 5 April 2020; Accepted: 6 May 2020; Published: 11 May 2020

**Abstract:** It has been hypothesized that solid tumors with invasive type of growth should possess intrinsic resistance to antiangiogenic therapy, which is aimed at cessation of the formation of new blood vessels and subsequent shortage of nutrient inflow to the tumor. In order to investigate this effect, a continuous mathematical model of tumor growth is developed, which considers variables of tumor cells, necrotic tissue, capillaries, and glucose as the crucial nutrient. The model accounts for the intrinsic motility of tumor cells and for the convective motion, arising due to their proliferation, thus allowing considering two types of tumor growth—invasive and compact—as well as their combination. Analytical estimations of tumor growth speed are obtained for compact and invasive tumors. They suggest that antiangiogenic therapy may provide a several times decrease of compact tumor growth speed, but the decrease of growth speed for invasive tumors should be only modest. These estimations are confirmed by numerical simulations, which further allow evaluating the effect of antiangiogenic therapy on tumors with mixed growth type and highlight the non-additive character of the two types of growth.

**Keywords:** mathematical oncology; spatially distributed modeling; reaction-diffusion-convection equations; computer experiment

**MSC:** 35K57; 35Q92; 92C05

#### **1. Introduction**

The use of mathematical methods has currently become a necessity in oncology. The identification of tumor cells [1], the design of nanomedical systems [2], and real-time adaptation of radiotherapy [3] are just a few examples of problems for which the solution already benefits from the use of simulation studies. One specific mathematical approach is the modeling of tumor growth and treatment, wherein a whole tumor and its microenvironment are considered as a single complex system. Its main goals are gaining new insights into various aspects of cancerous tumors [4,5] and the suggestions for the optimization of treatment protocols [6,7].

Clearly, a mathematical model of tumor growth must capture at least some of the most essential features of cancer cells' behavior and their interaction with the microenvironment. The most notable property of cancer cells is their ability for unlimited growth under favorable conditions [8]. However, tumor growth in tissue is restrained, first of all by the limited availability of nutrients. These aspects were taken into consideration already in the first non-spatially distributed phenomenological models of tumor growth [9,10]. Accounting for another hallmark of cancer—tissue invasion and metastasis—is possible in models that explicitly consider the spatial distribution of cancer cells. In continuous models, the intergrowth of tumor in normal tissue is usually realized via the introduction of parabolic terms, which include the intrinsic motility of tumor cells as a parameter [11,12]. As a rule, these models neglect the convective motion, which arises due to the proliferation of tumor cells. This process by itself can provide an increase in tumor volume even under zero cell motility. Accounting for it can be realized via hyperbolic equations and by itself results in a compact type of tumor growth [13,14]. Such a mechanism is often crucial for low-grade tumors, while the invasive type of growth usually begins to play an increasing role during their progression.

One more hallmark of cancer is sustained angiogenesis, i.e., the formation of new blood vessels that leads to an increase in tumor nutrient supply and thus promotes its growth. This process was incorporated in various ways in mathematical models of different complexity, from the ones governed by ordinary differential equations [15,16], to the complicated multiscale hybrid models [17,18].

The first antiangiogenic drug, bevacizumab, was approved for medical use in 2004 and is utilized currently along with about a dozen other angiogenesis inhibitors. Today, the majority of approved administration schemes, which include antiangiogenic drugs, combine them with various chemotherapeutic agents, which are aimed at direct cell killing [19]. Early experiments on mouse tumor models showed promising results regarding the use of antiangiogenic drugs in monotherapy regimes as well, since their administration as single agents allowed achieving significant delays in tumor growth. However, in the majority of the clinical trials, the administration of mono-AAT did not lead to any notable increase in survival [20]. Such a discrepancy is supposed to be linked to one obvious qualitative difference between preclinical and clinical tests: while the former were mostly conducted on localized primary tumors, the latter were mainly focused on the late-stage diseases [21]. Only relatively recently have clinical trials for early-stage localized tumors been initiated. Their preliminary data suggest the efficiency of mono-AAT for such tumors, at least in terms of tumor mass reduction [22].

Various mechanisms of resistance to treatment have been proposed in order to explain this effect [23]. One such mechanism is the accentuated invasiveness of tumor cells, which allows them to move away from nutrient-deprived zones. There exists a significant amount of experimental evidence that AAT accelerates tumor progression towards increasingly invasive phenotypes [21,24]. Based on this observation, it has been hypothesized that the tumors, which initially have an invasive phenotype, should possess intrinsic resistance to AAT, while compactly growing tumors should be more susceptible to it [23].

In this work, a continuous model of a monoclonal tumor growth in tissue is presented, which is able to reproduce this effect. The model allows accounting for both types of tumor growth, as well as their combination. Analytical estimations of tumor growth speed are obtained for tumors with both of the pure types of growth, which allow assessing the effect of antiangiogenic therapy on them. Numerical simulations are presented, which show good agreement with analytical estimations. They further allow evaluating the effect of antiangiogenic therapy on tumors with a mixed type of growth and provide insights into the mechanism of the interaction between the two types of growth, highlighting their non-additive character.

#### **2. Model**

#### *2.1. Equations*

The mathematical model of tumor growth, considered herein, represented a simplification of the model, previously developed by our research group. Its different versions were used for the investigation of several aspects of tumor growth and treatment [25–27]. There were four variables in this version of the model, which were a function of space and time coordinates, *x* and *t*: the density of tumor cells *n*(*x*, *t*), the fraction of necrotic tissue *m*(*x*, *t*), the concentration of glucose *g*(*x*, *t*), and the density of capillaries surface area *c*(*x*, *t*). All the variables, as well as all the used parameters should be strictly non-negative due to their physical meaning. The one-dimensional planar case was considered in this work, which was suitable for the consideration of a large spherically-symmetrical tumor. The following set of equations governs the dynamics of the model variables:

$$\begin{aligned} \text{tumor cells:} \quad \underbrace{\frac{\partial n}{\partial t}}\_{\partial t} &= \overbrace{Rn \cdot [1 - \sigma(g)] - Mn \cdot \sigma(g)}^{\text{prediction}} + D\_n \frac{\partial^2 n}{\partial x^2} - \frac{\partial (\text{In})}{\partial x}; \\ \text{increctic tissue:} \quad \underbrace{\frac{\partial m}{\partial t}}\_{\partial t} &= \overbrace{Mn \cdot \sigma'(g)}^{\text{death}} - \frac{\partial (\text{In})}{\partial x}; \\ \text{glucose:} \quad \underbrace{\frac{\partial g}{\partial t}}\_{\partial t} &= \overbrace{\frac{\partial m}{\partial t} - \frac{\text{cross}\_{\text{compensation}}}{\partial x \partial t}}^{\text{diffusion}} \underbrace{\text{diffusion}}\_{\left\{ \begin{aligned} \text{diffusion} \\ \text{degradation} \\ \text{degradation} \end{aligned} \right.} \end{aligned} \tag{1}$$

$$\text{explains:} \quad \underbrace{\frac{\partial c}{\partial t}}\_{\partial t} = \overbrace{-R[n+m]c};$$

$$\text{where} \quad \sigma(g) = \frac{1}{2} [1 - \tanh(\epsilon \{g - g\_{\mathcal{F}}\})];$$

$$\nabla \mathbf{I} = Bn \cdot [1 - \sigma(g)] + D\_n \frac{\partial^2 n}{\partial x^2}.$$

#### 2.1.1. Tumor Cells

The term tumor cell proliferation implies that it happens ceaselessly under a sufficient level of glucose, which was chosen to be the key nutrient, since it is indispensable for cell division [28]. With the fall of glucose level, the rate of proliferation slows down, and tumor cells die by necrosis. The drainage of necrotic tissue was neglected. The exact rates of the processes of cell proliferation and death were governed by a sigmoid function *σ*(*g*). Tumor cells were able to migrate throughout the tissue, which was governed by a diffusion-like term. The convective terms described the bulk motion of tissue elements, the velocity field **I** being determined by the dynamics of tumor cells. The expression for it was derived under the assumption of the constancy of the total density of tumor cells, necrotic tissue, and normal cells, which was normalized to unity. The normal cells were not considered in the model explicitly; however, it was assumed that the passive motion due to the arising convective flow was the only part of their dynamics.

#### 2.1.2. Glucose and Capillaries

The dynamics of glucose is comprised of its inflow from the capillaries into the tissue, consumption by tumor cells, and diffusion throughout the tissue. The inflow of glucose is governed primarily by the process of passive diffusion through the walls of capillaries [29]. Therefore, the rate of glucose inflow is proportional to the density of the capillaries' surface area and to the difference in glucose concentrations in blood and in tissue. Glucose concentration in blood was considered to be constant and was normalized to unity.

The capillary network degrades inside the tumor. This process has various reasons of a mechanical [30] and chemical nature [31], the details of which are difficult to account for in a reliable manner. Therefore, the degradation of capillaries was described by a rather phenomenological term. The volume of capillaries was considered to be negligible compared to the volume of cells, and therefore, their dynamics did not affect the convective velocity field. The convective motion of capillaries was also neglected. The normal density of capillaries' surface area was normalized to unity.

#### 2.1.3. Angiogenesis and Antiangiogenic Therapy

A more or less straightforward consideration of angiogenesis would require the introduction of an additional variable for the concentration of the main pro-angiogenic factor VEGF, which would affect the dynamics of capillaries. Such an approach has been utilized previously in different continuous models of tumor growth [12,32,33], including our models [27,34]. Herein, a more concise method was used for the sake of analytical study. It was assumed that during the growth of an untreated tumor, the concentration of VEGF was so high that it provided the maximum possible angiogenic effect throughout all the considered part of the tissue (the rate of this effect was limited by the number of receptors on endothelial cells). In fact, there are two separate effects, induced by VEGF, that affect the inflow of glucose from the capillaries and therefore needed to be considered in this model. One of them is the increase of the capillary density due to the formation of new capillaries. Another effect is the increase of the vascular permeability [35,36]. Due to the above-described assumption of the uniform action of VEGF, accounting for the latter effect came down to the increase in the value of parameter *P*, since it corresponded to the uniform identical increase in the permeability of all capillaries in the tissue. Assuming further that the local increase of the capillaries' surface area density, *c*, was proportional to its own local value throughout the tissue, one may conclude that for the consideration of the glucose inflow in tissue, an increase of *c* was equivalent to an analogous increase in the value of parameter *P*. Thus, the total effect of angiogenesis was accounted for in this model by the increase in the value of a single parameter *P*. Since AAT leads to the normalization of the structure of capillaries and further normalization of their density [37,38], its action as reflected by a decrease of *P* to a value, which was assumed to correspond to normal tissue.

#### *2.2. Parameters*

The parameters of the model were estimated according to the data of various experimental works. The basic set of parameters is listed in Table 1. The dimensionless model values of the parameters were the approximations of their normalized values, which were obtained with the use of the following normalization parameters: *tn* = 1 h for time, *xn* = 10−2 cm for length, *gn* = 1 mg/mL for glucose concentration, *cn* = 100 cm2/cm3 for normal capillary surface area density, which was close to its average value for human muscle [29], and *nn* = 3 × 108 cells/mL for the maximum density of tumor cells. The latter value was taken from the experimental work on the in vitro growth of multicellular tumor spheroids [39]. The values for the proliferation rate of tumor cells and their glucose consumption rate were also estimated according to the data of this work; however, it was assumed that these values should be proportionally diminished during the growth of a relevant tumor in tissue. The basic coefficient of tumor cells' motility corresponded to highly motile glioma cells [40]. This parameter is set to zero in Section 3.1 in order to focus on the compact type of tumor growth and is varied in Section 3.3 in order to consider tumors with a mixed type of growth. The value for capillaries' degradation rate was chosen in order to correspond to experimental observations, which showed that capillaries with adequate blood filling were very scarce inside the core of the tumors with radii of several millimeters [41]. The rate of death of tumor cells was taken to be significantly higher that that of their proliferation. The values of the parameters of function *σ*(*g*) could not be assessed straightforwardly, since the form of this function was chosen for phenomenological reasons. Therefore, it was merely taken to be rather close to a stepwise function, with the transition of tumor cells from proliferation to death happening at the concentrations of glucose an order of magnitude lower than that in blood. Importantly, the use of these values of variables *M*, *gcr*, and allowed keeping the concentrations of glucose positive within the necrotic zone.

**Table 1.** Basic set of model parameters.

In order to consider the effect of angiogenesis on tumor growth, the value of angiogenesis parameter *P* was varied up to a value ten times greater than the basic one. Such a limit was selected as a very approximate product of two estimated factors. Firstly, experiments in various mouse tumor models showed that the local density of microvessels near a tumor could increase three to six times [43]. Secondly, in the work [34], it was shown that a 2.5-fold increase in the permeability of tumor capillaries to glucose due to the action of VEGF was a physiologically reasonable estimation.

#### *2.3. Numerical Solving*

During numerical simulations, the set of equations Equation (1) was solved in a region with a size of several centimeters. The exact size *X* was adjusted for each case in order to be sufficiently small to spare computational time without imposing noticeable edge effects, while for all the variables, the zero-flux boundary condition was used at both edges. The convective flow speed was set to zero at the left boundary, where *x* = 0, resulting in the following equation for it:

$$\mathbf{I}(\mathbf{x},t) = \int\_0^\mathbf{x} Bn(r,t) \cdot [1 - \sigma(\mathcal{g}(r,t))] \mathrm{d}r + D\_{ll} \frac{\partial n(\mathbf{x},t)}{\partial \mathbf{x}}.$$

The following initial conditions were used, which represented a normal tissue with a small, 0.01 mm in width, colony of tumor cells, located near the left boundary, where capillaries were absent:

$$\begin{cases} n = 1, \\ m = 0, \\ \text{g} = 1, \\ c = 0 \end{cases} \quad \text{for } \mathbf{x} < = 0.1; \quad \begin{cases} n = 0, \\ m = 0, \\ \text{g} = 1, \\ c = 1 \end{cases} \quad \text{for } \mathbf{x} > 0.1. \tag{2}$$

The method of splitting into physical processes was used for all variables, i.e., kinetic equations, diffusion equations, and convective equations were solved successively during each time step. The implicit Crank–Nicholson scheme was used for the glucose diffusion equation. Since the glucose diffusion term provided maximum rate of local change among all the variables and thus required a sufficiently small time step even for implicit solving, it was decided to solve the cell migration equation by a simpler explicit forward Euler scheme, and all kinetic equations by the explicit Euler method. Convective equations were solved using the flux-corrected transport algorithm with explicit anti-diffusion stage. The last method was introduced in the work [44], while other classical methods were described in many books (see, e.g., [45]). The choice of time and space steps is justified in the following sections. For optimization purposes, the function *σ*(*g*) was not recalculated every time; instead, it was precalculated for about ten thousand values of *g*, evenly distributed on the segment [0, 1], and only thusly obtained values were used during the calculations as approximations of the actual values of *σ*(*g*).

The tumor growth speed *Vgr*(*t*) was calculated as the rate of change of tumor radius, which was evaluated as the maximum space coordinate, at which *n* ≥ 0.1. In all the simulations, after an initial transitional period, the speed of tumor growth seemed to tend asymptotically to some constant value. In order to ensure that this value could be estimated with suitable precision, for each simulation, the part of the function *V*(*t*), obtained after the manually designated initial transitional period, was fitted via the least-squares method by a function of the form:

$$V\_{\mathbb{S}^r}(t) = V + \sum\_{i=1}^k a\_i e^{-b\_i t} \mu\_i$$

wherein the value *k* = 3 was selected manually as the one providing sufficiently fine approximation. A simulation was stopped manually if the current tumor growth speed and *V*, expressed in mm/week, were equal up to three decimal places. If it did not happen until the proximity of the tumor to the right boundary notably affected its growth speed, the simulation was rerun in a larger region. Further throughout the text, this limiting value *V* is meant for tumor growth speed.

For some simulations, the number of tumor cells was estimated analogically as:

$$N = \lim\_{t \to \infty} \int\_0^X n(x) dx\_\prime$$

the integral being calculated numerically by the elementary rectangle rule.

The computational codes were implemented in C++ and can be found in the Supplementary Materials.

#### **3. Results**

#### *3.1. Compact Type of Growth*

At first, let us obtain an analytical estimation for the growth speed of a solid tumor with zero cell motility, i.e., *Dn* = 0, which thus has a purely compact type of growth. For this purpose, let us seek the solution of Equation (1) in the form of a wave, traveling with constant shape and speed *V* in an infinite region. In such a case, the governing equations can be reduced to a system of ordinary differential equations by introducing the traveling coordinate frame *z* = *x* − *Vt*:

$$\begin{aligned} \text{tumor cells:} \quad &Bu \cdot [1 - \sigma(\mathfrak{g})] - Mn \cdot \sigma(\mathfrak{g}) - \frac{\partial(\mathrm{Im})}{\partial z} + V \frac{\partial n}{\partial z} = 0;\\ \text{reccctic tissue:} \quad Μ \cdot \sigma(\mathfrak{g}) - \frac{\partial(\mathrm{Im})}{\partial z} + V \frac{\partial m}{\partial z} = 0;\\ \text{glucose:} \quad &Pc[1 - \mathfrak{g}] - Qu \cdot [1 - \sigma(\mathfrak{g})] + D\_{\mathbb{F}} \frac{\partial^2 \mathfrak{g}}{\partial z^2} + V \frac{\partial \mathfrak{g}}{\partial z} = 0;\\ \text{capillaries:} \quad & -R[n + m]c + V \frac{\partial c}{\partial z} = 0; \end{aligned}$$

$$\begin{array}{ll}\text{where} & \sigma(\mathcal{g}) = \frac{1}{2} [1 - \tanh(\epsilon \{\mathcal{g} - \mathcal{g}\_{\mathcal{F}}\})]; \\ & \mathbf{I}(z) = \int\_{-\infty}^{z} Bn(r) \cdot [1 - \sigma(\mathcal{g}(r))] \mathbf{d}r; \\ & V = \lim\_{z \to +\infty} \mathbf{I}(z). \end{array}$$

The living part of the tumor was considered as a planar front, which propagated towards the right boundary, which represented the normal tissue, and left the necrotic zone behind it. That is formalized by the following boundary conditions:

$$\begin{cases} \begin{array}{l} n=0, \\ m=1, \\ \frac{\partial \chi}{\partial z} = 0, \\ c=0, \\ \mathbf{I}=0 \end{array} \end{cases} \quad \begin{cases} n=0, \\ m=0, \\ \mathbf{g}=1, \\ c=1, \\ \frac{\partial \mathbf{I}}{\partial \overline{z}} = 0, \end{cases} \quad \text{for } z \to +\infty. \tag{4}$$

Further, in order to obtain an analytically tractable case, let us consider the following limitations:

The last condition is sensible only if the tumor has a clear boundary, i.e.,

$$\exists \mathfrak{x}\_b \colon \, \forall \mathfrak{x} > \mathfrak{x}\_b \colon n(\mathfrak{x}) + m(\mathfrak{x}) = 0,$$

These limitations result in the following system:

$$\begin{aligned} \text{tumor cells:} \quad & \begin{cases} n = 0 \text{ if } g \le g\_{cr}, \\ Vn' - (\text{In})' + Bn = 0 \text{ if } g > g\_{cr} \end{cases} \\ \text{necrotic tissue:} \quad & m = \begin{cases} 1 \text{ if } g \le g\_{cr} \\ 0 \text{ if } g > g\_{cr} \end{cases} \\ \text{glucose:} \quad & D\_{\mathcal{S}} \mathcal{Y}^{\prime\prime} + V \mathcal{Y}^{\prime} + \text{Pc} [1 - \mathcal{g}] - \mathcal{Q}n = 0; \\ \text{capillaries:} \quad & c = \begin{cases} 1 \text{ if } n + m = 0, \\ 0 \text{ if } n + m > 0; \end{cases} \end{aligned} \tag{5}$$

$$where \quad \mathbf{I}(z) = \int\_{-\infty}^{z} Bn(r) \text{d}r; \quad V = \lim\_{z \to +\infty} \mathbf{I}(z),$$

where primes denote differentiation with respect to *z*.

Finally, let us specify the form of the tumor cells' distribution in the sought solution as a limiting case of a piecewise function, which can be equal only to zero and one, like the functions for necrotic tissue and capillaries' distribution. Thus, the sought solution has the form depicted in Figure 1, where the origin of the z-axis is placed at the front of the tumor for convenience. This solution can be split up into three regions: (1) necrotic core, (2) proliferating rim of yet unknown width *L*, in which *n* = 1 and *g* > *gcr*, and (3) normal tissue. The expression for the tumor growth speed is now simplified to:

−∞

$$V = \mathcal{I}(0) = BL\dots$$

**Figure 1.** The form of the solution of Equation (5), which is searched for analytically.

The distribution of glucose has to be found in order to estimate tumor growth speed. This procedure, described in Appendix A, yields the following implicit expression for *L*:

$$[1 - \frac{\mathcal{Q}}{B} + \frac{\mathcal{Q}D\_{\mathcal{S}}}{B^2 L^2} [1 - e^{-BL^2/D\_{\mathcal{S}}}] [1 - \frac{2}{1 + \sqrt{1 + 4D\_{\mathcal{S}} P/(B^2 L^2)}}] = \mathcal{g}\_{cr} \tag{6}$$

Its numerical solving resulted in *L* ≈ 1.325 under the values of the parameters used, taken from the basic parameter set (see Table 1), which led to the tumor growth speed *V* ≈ 0.223 mm/week. There was under a ten-fold greater value for angiogenesis parameter, i.e., *P* = 40, *L* ≈ 2.602 and *V* ≈ 0.437 mm/week, which meant almost a two-fold increase in tumor growth speed. These values were in good correspondence with the speeds of the in vitro growth of multicellular tumor spheroids, experimentally obtained in the work [39] (≈0.5 mm/week), based on which several values of the model parameters were estimated.

The estimation of *V* can be significantly simplified due to the presence of a small parameter:

$$\delta = \frac{BL^2}{D\_{\S}} = \frac{V}{D\_{\S}/L}\prime$$

the smallness of which is due to the fact that the tumor growth speed should be significantly lower than the characteristic speed of glucose diffusion within the proliferating rim. For example, for the values of parameters from their basic set, *δ* ≈ 1.7 × 10−4. Expansion of the expression (6) in a Taylor series up to *o*(*δ*) yields:

$$\mathbf{g}\_{\sigma\tau} = 1 - \frac{\mathcal{Q}}{B} + \frac{\mathcal{Q}}{B} [1 - \sqrt{\frac{B}{P}} \sqrt{\delta} - \frac{\delta}{2} (1 - \frac{B}{P}) + o(\delta)].$$

Moreover, it is convenient to neglect *B*/*P* as a small parameter compared to one, e.g., *B*/*P* ≈ 2.5 × 10−3 for the basic values of the model parameters. This allows obtaining the equation for an approximate value of proliferating rim width, *L*:

$$g\_{cr} = 1 - \frac{Q}{\sqrt{PD\_{\mathcal{S}}}} \tilde{L} - \frac{Q}{2D\_{\mathcal{S}}} \tilde{L}^2 \tilde{\lambda}$$

which can be solved, only its positive root having physical meaning. That leads to the following formula for approximate tumor growth speed:

$$\vec{V} = B\vec{L} = B\sqrt{D\_{\mathcal{S}}} \left\{ \sqrt{\frac{1}{P} + 2\frac{1 - \mathcal{g}\_{cr}}{Q}} - \sqrt{\frac{1}{P}} \right\},\tag{7}$$

which under the considered parameter values provides the values of tumor growth speed, equal up to three decimal places to the ones derived numerically from Equation (6). Note that this expression has a limit under *P* → ∞:

$$
\tilde{\nu}\_{lim} = B \sqrt{\frac{2D\_{\tilde{\mathcal{S}}}}{Q} (1 - g\_{cr})} \tilde{\nu}
$$

that is ≈0.651 mm/week under the basic set of parameters. Thus, the speed of growth of a considered tumor cannot increase more than three-fold due to angiogenesis. Under *P* = 40, it is about two-thirds of its limit value, while under *P* = 80, it would be around three-quarters of it.

These estimations allow suggesting that for compactly growing tumors, AAT may provide a several times decrease of their growth speed. Of note, there is no account for the drainage of necrotic tissue in the model, while its consideration should enhance the result. Moreover, in this case, AAT might be able to lead to a complete tumor growth stop, as well as its shrinkage, which is sometimes observed experimentally [46].

The estimations of tumor growth speed under finite values of , *M*, and *R* are quite difficult if at all impossible to perform analytically, and they were performed via numerical simulations. However, at first, it was examined how well the results, obtained in a numerical experiment, could correspond to the already obtained analytical estimations. For this purpose, the sets of four simulations with different time and space steps were conducted for each of the six values of the angiogenesis parameter within the considered range *P* ∈ [4, 40]. In each of the simulations, the time step *τ* and space step *h*

related as *τ* = *h*2. The values of *R* and *M* were chosen to be close to the maximum values allowed by numerical calculations, i.e., *R* = *M* = 0.1/*τ*. The value of was chosen to be 105, which practically resulted in a stepwise function *σ*(*g*) due to the method of its implementation in the code.

The dots in Figure 2a denote the values of tumor growth speed, obtained numerically in these simulations. The analytically obtained values are designated by crosses of the corresponding colors on the vertical axis. For every used value of *P*, the best fit quadratic function was built in order to approximate the numerical value for the tumor growth speed under *h* → 0, *τ* → 0, *M* → ∞, and *R* → ∞. The thusly obtained values were in a good correspondence with the analytical ones, being only slightly smaller than them. The discrepancy between the values increased with the increase of *P*, constituting less than 1% under *P* = 4 and less than 5% under *P* = 40. Figure 2b–e shows the distributions of the model variables on the 20th day of simulations under the values of parameters, designated by the corresponding letters in Figure 2a. They showed that, expectedly, the numerically obtained profiles of tumor cells differed from stepwise functions; however, their front edges were getting steeper under finer discretization. Of note, the tumors that were obtained under coarser discretization, shown in Figure 2c,e, had greater radii on the 20th day than the corresponding tumors obtained under finer discretization, shown in Figure 2b,d. This did not reflect the relation between the growth speeds of the tumors in these simulations, since these profiles were obtained during the transitional period of tumor growth. As Figure 2a shows, tumor growth speed changed differently with the refinement of the discretization under different values of *P*.

**Figure 2.** (**a**) The dots denote the values of tumor growth speed *V*, obtained numerically in the simulations of Equation (1) under designated values of angiogenesis parameter *P* and space step *h*, with time step *τ* = *h*2, *Dn* = 0, *R* = *M* = 0.1/*τ*, = 105 and the values of other parameters taken from the basic set. The dashed lines represent the best fit quadratic functions for every value of *P*. The crosses on the vertical axis mark the analytically derived values of tumor growth speed for the values of *P* used. (**b**–**e**) The distributions of the model variables on the 20th day of simulations under the values of the parameters, designated by the corresponding letters in (**a**).

The analogical results, obtained in simulations with the values of parameters , *R*, and *M*, taken from the basic set, are shown in Figure 3. With the use of the basic values for these parameter, the approximated numerical values for tumor growth speed increased modestly for each considered value of *P*, becoming ≈1.8–2.4% higher than the corresponding analytically estimated values. Certainly, a further decrease of either or *R*, as well as the increase of either proliferation rate *B*, or glucose consumption rate *Q*, or glucose diffusion coefficient *Dg*, would lead to a further increase in tumor growth speed. However, too large alterations of the values of these parameters would lead to physically meaningless results involving negative glucose concentrations and/or explosive tumor growth, thus indicating the limitations of the presented model. Of note, the decrease of *M* by itself had only a small effect on the tumor growth speed, as well as on the profile of glucose concentration, but its change largely influenced the border between the profiles of necrotic tissue and tumor cells, which were considered to be yet alive.

**Figure 3.** (**a**) The dots denote the values of tumor growth speed *V*, obtained numerically in the simulations of Equation (1) under designated values of angiogenesis parameter *P* and space step *h*, with time step *τ* = *h*2, *Dn* = 0, and the values of other parameters taken from the basic set. The dashed lines represent the best fit quadratic functions for every value of *P*. The crosses on the vertical axis mark the analytically derived values of tumor growth speed for the values of *P* used. (**b**–**e**) The distributions of the model variables on the 20th day of simulations under the values of the parameters, designated by the corresponding letters in (**a**).

#### *3.2. Invasive Type of Growth*

Now, let us obtain an analytical estimation for the growth speed of a solid tumor with non-zero cell motility *Dn*, but neglecting the bulk motion of tissue elements, expressed in Equation (1) by hyperbolic terms. Once again, the solution was sought in the form of a wave, traveling with a constant shape and speed *V* in an infinite region. The boundary conditions are now modified as:

$$\begin{cases} \begin{array}{l} n=0, \\ m=m^\*, \\ \text{g}=\text{g}^\*, \\ c=0 \end{array} \end{cases} \quad for \; z \to -\infty; \quad \begin{cases} \begin{array}{l} n=0, \\ m=0, \\ \text{g}=1, \\ c=1 \end{array} \quad for \; z \to +\infty, \tag{8}$$

where *g*∗ and *m*∗ are introduced as the limiting constant values of glucose concentration and necrotic tissue fraction at *z* → −∞. For the estimation of *V*, it is necessary to consider the asymptotic behavior of the solutions at *z* = ±∞. For this purpose, let us linearize the system (1) at the boundary values (8), neglecting the convective terms, and thus obtain two systems of ODEs with constant coefficients:

$$\begin{cases} D\_n n\_-^{\prime\prime} + V n\_-^{\prime} + B n\_- \cdot [1 - \sigma(\mathfrak{g}^\*)] - M n\_- \cdot \sigma(\mathfrak{g}^\*) = 0; \\ V m\_-^{\prime} + M n\_- \cdot \sigma(\mathfrak{g}^\*) = 0; \\ D\_\mathfrak{g} \mathfrak{g}\_-^{\prime\prime} + V \mathfrak{g}\_-^{\prime} + P c\_- [1 - \mathfrak{g}^\*] - Q n\_- \cdot [1 - \sigma(\mathfrak{g}^\*)] = 0; \\ V c\_-^{\prime} - R m^\* c\_- = 0. \end{cases}$$

and:

$$\begin{cases} D\_n n\_+^\prime + V n\_+^\prime + B n\_+ \cdot [1 - \sigma(1)] - M n\_+ \cdot \sigma(1) = 0; \\ V m\_+^\prime + M n\_+ \cdot \sigma(1) = 0; \\ D\_\mathcal{S} \mathcal{g}\_+^\prime + V \mathcal{g}\_+^\prime - P \mathcal{g}\_+ - Q n\_+ \cdot [1 - \sigma(1)] = 0; \\ V c\_+^\prime - R [n\_+ + m\_+] = 0, \end{cases}$$

where primes denote differentiation with respect to *z* = *x* − *Vt* and subscripts − and + refer to the linearized problems at *z* → −∞ and *z* → +∞, respectively. The solutions to these problems should be sought in the form (*n* ±, *n*±, *m*±, *g* ±, *g*±, *c*±)*T* ∼ **k**±*exp*(*μ*±*z*), which reduces these systems of linear differential equations to eigenvalue problems.

For *z* → −∞, the eigenvalues are:

$$\begin{array}{c} \mu\_{1,2}^{-} = 0, \ \mu\_{3}^{-} = Rm^{-}/V, \ \mu\_{4}^{-} = -V/D\_{\mathbb{S}'}\\ \mu\_{5,6}^{-} = [-V/2 \pm \sqrt{V^{2}/4 + D\_{\text{fl}}\{M\sigma(\text{g}^{-}) - B[1 - \sigma(\text{g}^{-})] \}}]/D\_{\text{th}} \end{array}$$

where *μ*− 3 is positive, *μ*− 4 is negative, and *μ*− 5,6 have opposite signs, since *Mσ*(*g*−) − *B*[1 − *σ*(*g*−)] > 0, i.e., the death rate exceeds the proliferation rate at sufficiently low values of *z*, which follows from the physical meaning of the model.

For *z* → +∞, the eigenvalues are:

$$\begin{aligned} \mu\_{1,2}^{+} &= 0, \; \mu\_{3,4}^{+} = \left[ -V/2 \pm \sqrt{V^2/4 + D\_{\mathcal{S}}P} \right] / D\_{\mathcal{S}'}\\ \mu\_{5,6}^{+} &= \left[ -V/2 \pm \sqrt{V^2/4 - D\_{\mathcal{n}}\{B[1 - \sigma(1)] - M\sigma(1)\}} \right] / D\_{\mathcal{n}'} \end{aligned}$$

where *μ*+ 3,4 have opposite signs and, depending on the value of *V*, *μ*+ 5,6 are either complex numbers with negative real parts, or real numbers with opposite signs, or both are equal to zero. Therefore, for sufficiently low values of *V*, the solution may oscillate for large values of *z*, yielding regions with negative values of *n* and *m*, which is physically unrealistic. Thus, the following restriction on tumor growth speed is sufficient for physically reasonable solutions:

$$V \ge V\_{\min} = 2\sqrt{D\_n\{B[1-\sigma(1)] - M\sigma(1)\}} \approx 2\sqrt{B D\_n} \tag{9}$$

the last two values being extremely close under the basic set of parameters, differing by less than 10−13% even under a decrease of up to the value of 20. For the basic set of parameters, *Vmin* ≈ 1.063 mm/week, which was in a good correspondence with the speed of growth of highly invasive tumors [40]. Formula (9) corresponds well to the well-known result regarding Fisher's equation, which can be written out in the notation used herein as:

$$\frac{\partial n}{\partial t} = Bn(1 - n) + D\_n \frac{\partial^2 n}{\partial x^2}. \tag{10}$$

For this equation, the range of the speed of monotone waves satisfies:

$$V \ge V\_{\text{min}}^f = 2\sqrt{B D\_{\text{ref}}}$$

this result not being affected by the alteration of the proliferation term as long as the proliferation rate tends to *B* for *n* → 0. For Fisher's equation, it is known that if its initial conditions have a compact support, i.e.,

$$n(\mathbf{x},0) = n\_0(\mathbf{x}) \ge 0, \ n\_0(\mathbf{x}) = \begin{cases} 1 & \text{if } \mathbf{x} \le \mathbf{x}\_1 \\ 0 & \text{if } \mathbf{x} \ge \mathbf{x}\_2 \end{cases} \tag{11}$$

where *x*1 < *x*2 and *n*0(*x*) is continuous in *x*1 < *x* < *x*2, then the solution *n*(*x*, *t*) of Equation (10) evolves to a traveling wavefront solution with the speed *V f min* = 2 √*BDn* [47].

The numerical simulations of the system (1) under the neglect of convective terms with the initial conditions (2), which represent a function with a compact support, show that they evolve to the wavefronts, traveling with speeds very close to 2√*BDn*. For the basic set of parameters, numerical simulations gave *V* ≈ 1.061 mm/week, yielding the same tumor growth speed up to three decimal places (as well as the same number of cells *N* ≈ 1.383) under space steps *h* = {0.2, 0.1, 0.05, 0.025} and time steps *τ* = *h*2. Thus, much coarser discretization could provide more accurate results for the considered system under the neglect of convective terms. A dozen simulations under other values of parameters, chosen randomly from the ranges *B* ∈ [0.05, 0.3] and *Dn* ∈ [0.001, 0.4], as well produced the values of tumor growth speeds that differed from 2√*BDn* by no more than 0.5%.

Since the considered system could be viewed as an augmentation of Fisher's model, the obtained results suggested that the introduced modifications should not affect the wave speed of its solution, which evolved from a function having a compact support, or at least they should lead to only a very small correction of it within the physiologically justified range of parameters. Most importantly, the parameter *P* is not present in Formula (9), which implies that angiogenesis should not influence the growth speed of a considered invasive tumor at all. Numerical modeling speaks in favor of this result. The simulations under six values of the angiogenesis parameter within the considered range *P* ∈ [4, 40] yielded tumor growth speeds equal up to three decimal places, which is illustrated by Figure 4a. The simulations were run with space step *h* = 0.1 and time step *τ* = *h*2 = 0.01, which nevertheless did not allow saving a sufficient amount of computational time, compared to the simulations of Section 3.1, since in this case, the tumor growth speed tended sufficiently slower to a constant value, thus requiring longer simulation on larger domains.

Of note, the number of tumor cells notably grew with the increase of *P*, as Figure 4a shows, resulting in an ≈72% increase under 10-fold magnification of the basic value of *P* = 4. This effect was well noticeable when comparing the profiles of tumor cells under *P* = 4 and *P* = 40, which are depicted in Figure 4b,c on the 20th day of tumors growth.

**Figure 4.** (**a**) The red dots denote the values of tumor growth speed *V*; the green dots denote the values of number of tumor cells *N*, which were obtained numerically in the simulations of Equation (1) with the neglect of hyperbolic terms (i.e., **I** = 0) under designated values of angiogenesis parameter *P*, values of other parameters, taken from the basic set, space step *h* = 0.1, and time step *τ* = 0.01. The dashed lines are interpolations of data points. (**b**,**c**) The distributions of the model variables on the 20th day of simulations under *P* = 40 and *P* = 4 correspondingly.

#### *3.3. Mixed Type of Growth*

The conclusion of the indifference of invasive tumor growth speed to angiogenesis was valid under the full neglect of the convective component of tumor growth. However, it was reasonable to assume that simultaneously accounting for it should nevertheless lead to a non-zero effect of AAT on the growth speed of an invasive tumor, since the convective component of tumor growth should be affected by the number of its cells. In order to evaluate the effect of antiangiogenic therapy on tumors with a mixed type of growth, the sets of six simulations under different values of angiogenesis parameter *P* were conducted for eight values of tumor cell motility *Dn*. As a compromise between computing resources and accuracy, the values of space step *h* = 0.01 and time step *τ* = 10−4 were used. The obtained tumor growth speeds are designated in Figure 5a, where the data for the case with zero cell motility, already shown in Figure 3a, is added.

Let introduce the parameter of "maximum antiangiogenic effect" as:

$$(V|\_{P=40} - V|\_{P=4}) / V|\_{P=40}$$

which indicates how the tumor growth speed would slow down upon a decrease of angiogenesis parameter *P* from 40 to four, the latter corresponding to the normal microvasculature level. Its dependence on tumor cell motility is shown in Figure 5b, while under *Dn* = 0, it was ≈48%. The graph suggests that the effect of mono-AAT should inversely correlate with the invasiveness of the tumor; however, it should be non-zero even for highly invasive tumors, since the maximum antiangiogenic effect, obtained under high values of tumor cell motility, *Dn* = 0.1 and *Dn* = 0.33, was ≈13% and ≈11%, correspondingly.

**Figure 5.** (**a**) The dots denote the values of tumor growth speeds *V*, obtained numerically in the simulations of Equation (1) with space step *h* = 0.01 and time step *τ* = 10−4, under designated values of angiogenesis parameter *P* and tumor cell motility *Dn* and the values of other parameters taken from the basic set. The dashed lines are interpolations of data points. (**b**) The graph of "maximum angiogenic effect", based on the data, shown in (**a**), which shows how the tumor growth speed would slow down upon the decrease of angiogenesis parameter *P* from 40 to four under designated values of *Dn*.

An interesting property of the considered tumors with a mixed type of growth was the fact that their growth speed was not equal to the sum of two speeds: *Vcomp* that was obtained in the absence of cell motility—i.e., during purely compact growth—and *Vinv* ≈ 2 √*BDn* that was obtained in absence of convective flows—i.e., during purely invasive growth. This property is illustrated in Figure 6, which shows the dependence of the parameter (*V* − *Vcomp*)/*Vinv* on *P* for different *Dn*. This parameter may be treated as the ratio of the increase in tumor growth speed, caused by the mobilization of initially immotile cells, to the speed of pure invasive growth itself. The reason for the non-additivity of the growth speeds was the fact that the redistribution of the tumor cells, caused by their migration, led to the change in the number of cells in the proliferating state and altered the convective velocity field. This effect had an ambiguous character. Firstly, the protrusion of tumor cells' profile towards the region with a normal density of microvasculature led to the increase in the effective rate of capillaries' degradation, which resulted in a decrease of glucose inflow to the tumor and a subsequent decrease of the pool of proliferating cells. Secondly, this protrusion itself enlarged the pool of tumor cells, located in the region with a concentration of glucose sufficient for their proliferation. The first aspect dominated under small values of *Dn*, leading to a much smaller increase in tumor growth speed, than 2 √*BDn*. This effect was more pronounced under high values of *P*, since in this case, tumors had more cells, and therefore enlarged sizes of such protrusions that led to accelerated degradation of capillaries. According to the simulations, a tumor with motile cells may even grow slower than the same tumor with immobilized cells. The second aspect dominated under high values of *Dn*, resulting in the fact that the tumor growth speed was higher than *Vconv* + *Vdi f* under *Dn* = 0.33 for all the considered values of *P*.

**Figure 6.** The dots denote the ratios of the differences of tumor growth speed *V* and *Vcomp*, which is the growth speed of the same tumor under zero cell motility, to *Vinv*, which is the growth speed of the same tumor under the neglect of convective flows. All values were obtained numerically in the simulations of Equation (1) with space step *h* = 0.01 and time step *τ* = 10−4 under designated values of angiogenesis parameter *P* and tumor cell motility *Dn* and the values of other parameters taken from the basic set. The dashed lines are interpolations of data points.

#### **4. Discussion**

The goal of this paper was to study, by means of mathematical modeling, the effect of mono-antiangiogenic therapy (AAT) on monoclonal tumors, which have different types of growth. The analytical estimations performed allowed suggesting that AAT may provide a several times decrease of growth speed for compact tumors, but the decrease of growth speed for invasive tumors should be much less significant. Thus, the tumors, which initially had an invasive phenotype, should possess intrinsic resistance to AAT, while compactly growing tumors should be more susceptible to it. This conclusion corresponds well to the preclinical and clinical data [21,23]. To the best of my knowledge, this effect was not reproduced before via a continuous model of tumor growth, while this ability of the model should be crucial for the investigation of AAT.

The presented model was built on the fundamental principle that the diffusion and consumption of nutrients limit the growth of a solid tumor. This principle underlines a multitude of existing continuous models that account for spatio-temporal interactions between tumor cells and nutrients (see, e.g., [48] for a review). The main distinguishing feature of the presented model was the simultaneous accounting of the intrinsic motility of malignant cells and the convective motion, arising due to their proliferation, which allowed considering two types of tumor growth, as well as their combination. Despite the relative simplicity that this approach has under the assumption of the spherical symmetry of a tumor, it is relatively rarely used, and usually one of the two types of tumor growth is ignored in continuous models. Some of the few exceptions were the work [49], which provided a theoretical justification for the phenomenon of the dominance of a metastatically active population in an avascular heterogeneous tumor, and the work [50], which explained from a theoretical point of view the effect of tumor cells' migration within a multicellular tumor spheroid. The works [51,52], which also utilized such an approach, were devoted to the investigation of the fingering instability of the avascular tumor surface in the two-dimensional case. This instability could be regarded as a way for a tumor to counteract the diffusional limitations of nutrient inflow. While this phenomenon may affect the tumor growth speed quantitatively, the main qualitative results of the presented work, obtained under the assumption of radial symmetry, could hardly be affected in higher dimensions.

The presented model accounts for glucose as the main nutrient, since it is an indispensable substrate for cell division [28], while most of the relevant models consider tumor cell proliferation to be dependent on oxygen [48]. While such a difference may not be of major importance to the model presented herein, as well as to the models of avascular tumor growth, it should be noted that the approach for modeling the dynamics of oxygen with implicit accounting of tumor vasculature should differ from that for glucose due to the specific features of their blood transport, transvascular transport, and tumor metabolism, which were discussed in the previous paper by our research group [34].

The account of angiogenesis, utilized herein, is rather schematic; however, it allows considering the two major effects of AAT on a qualitative level, allowing performing analytical estimations that suit the purposes of this work. Under more straightforward consideration of angiogenesis, i.e., the introduction of a separate variable for VEGF, which would affect the dynamics of capillaries, one would obtain a distribution of capillaries with non-uniform alterations in their density and permeability [27,34]. However, for every simulation of such an extended model, there will exist a corresponding value of *P*, which would provide the same increase in tumor growth speed in the simplified model presented herein. Interestingly, the analytical expression for the speed of growth of a compact tumor with no capillaries inside it allowed obtaining a limit of its growth speed under the infinite increase in the number of capillaries and/or in their permeability. However, the estimations with physiologically based values of parameters suggested that achieving more than ≈ 80% of this limit value of speed was unlikely for real tumors. The analytical expression for the speed of growth of invasive tumor indicated that angiogenesis did not influence it at all under the neglect of the convective component of tumor growth, which was confirmed by numerical simulations. However, the simulations, which accounted for both reasons of growth, suggested that in terms of tumor speed reduction, the maximum possible angiogenic effect for highly invasive monoclonal tumors should be around 10–15%.

Importantly, the neglect of convective motion would have led to the misleading conclusion of the indifference of tumor growth speed to angiogenesis for low-invasive tumors as well, since it was not affected by the exact value of tumor cell motility. It should be noted that the accounting of the convective component of tumor growth should be also crucial for modeling of other types of antitumor therapy, since the decrease in the number of tumor cells, caused by any treatment, in reaction-diffusion models should lead to an underestimated decrease of tumor growth speed.

The numerical simulations performed highlighted the non-additive character of the two types of tumor growth. Namely, the addition of a small enough motility to initially immobilized tumor cells should lead to a notably smaller increase in tumor growth speed that could be expected from the analytical estimations of the speed of tumor invasion. Interestingly, the simulations suggested that a low-invasive tumor may even grow slower than the same tumor with immobilized cells. Under sufficiently high cell motility, on the contrary, the two types of tumor growth produced a synergistic effect.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2227-7390/8/5/760/s1.

**Funding:** The publication has been prepared with the support of the "RUDN University Program 5-100".

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **Appendix A. Analytical Estimation of Compact Tumor Growth Speed**

The distribution of glucose in Figure 1 is searched as a piecewise function that must be continuously differentiable for the continuity of glucose concentration and its flow:

$$\begin{array}{c} \text{g}(z) = \text{g}\_{I}(z), \; z < -L; \\ \text{g}(z) = \text{g}\_{II}(z), \; -L < z < 0; \\ \text{g}(z) = \text{g}\_{III}(z), \; 0 < z; \end{array}$$

$$\begin{array}{c} \mathbf{g}(-L) = \mathbf{g}\_{I}(-L) = \mathbf{g}\_{II}(-L) = \mathbf{g}\_{cr}, \ \mathbf{g}\_{I}'(-L) = \mathbf{g}\_{II}'(-L);\\ \mathbf{g}(0) = \mathbf{g}\_{II}(0) = \mathbf{g}\_{III}(0), \ \mathbf{g}\_{II}'(0) = \mathbf{g}\_{III}'(0). \end{array}$$

Within the necrotic core, the equation for glucose turns into:

$$D\_{\mathfrak{K}} \mathfrak{g}\_I'' + BL\mathfrak{g}\_I' = 0\_{\prime\prime}$$

the general solution of which is:

$$\mathbf{g}\_I = \mathbf{C}\_1^I + \mathbf{C}\_2^I e^{-BLz/D\_{\mathcal{S}\_{\mathbf{s}'}}} $$

for which the boundary conditions allow only for:

$$\mathcal{g}\_I = \mathcal{g}\_{cr}.$$

Within the proliferating rim, the equation for glucose is converted into:

$$-Q + D\_{\mathcal{S}} \mathcal{g}\_{II}^{\prime\prime} + BL\_{\mathcal{S}} \mathcal{g}\_{II}^{\prime} = 0,$$

the general solution of which is:

$$\lg\_{II} = \mathbf{C}\_1^{II} + \mathbf{C}\_2^{II} \mathbf{e}^{-BLz/D\_{\mathcal{S}}} + \frac{Qz}{BL} \mathbf{J}$$

Stitching of *gI* and *gI I* at *z* = −*L* yields:

$$\begin{aligned} \mathbf{g}\_I(-L) &= \mathbf{g}\_{II}(-L): & \mathbf{g}\_{\mathcal{F}} &= \mathbf{C}\_1^{II} + \mathbf{C}\_2^{II} \mathbf{e}^{BL^2/D\_{\mathcal{S}}} - \frac{Q}{B}; \\ \mathbf{g}\_I'(-L) &= \mathbf{g}\_{II}'(-L): & 0 &= -\frac{BL}{D\_{\mathcal{S}}} \mathbf{C}\_2^{II} \mathbf{e}^{BL^2/D\_{\mathcal{S}}} + \frac{Q}{BL}. \end{aligned}$$

Therefore:

$$\mathbb{C}\_1^{II} = \mathbb{g}\_{cr} + \frac{Q}{B} [1 - \frac{D\_{\mathbb{g}}}{BL^2}]; \ \mathbb{C}\_2^{II} = \frac{QD\_{\mathbb{g}}}{B^2 L^2} e^{-BL^2/D\_{\mathbb{g}}}.$$

For normal tissue, the glucose equation transforms into:

$$P[1 - \mathcal{g}\_{III}] + D\_{\mathcal{S}} \mathcal{g}\_{III}^{\prime\prime} + BL\mathcal{g}\_{III}^{\prime} = 0,$$

the general solution of which is:

$$\mathbf{g}\_{III} = \mathbf{1} + \mathbf{C}\_1^{III} \exp(\frac{-BL - \sqrt{B^2 L^2 + 4D\_\mathcal{S}P}}{2D\_\mathcal{S}} z) + \mathbf{C}\_2^{III} \exp(\frac{-BL + \sqrt{B^2 L^2 + 4D\_\mathcal{S}P}}{2D\_\mathcal{S}} z).$$

Since the glucose concentration is limited at *z* → +∞, then *CIII* 2 = 0. Stitching of *gI I* and *gIII* at *z* = 0 results in:

$$\begin{aligned} g\_{II}(0) &= g\_{III}(0): & \mathcal{C}\_1^{II} + \mathcal{C}\_2^{II} &= 1 + \mathcal{C}\_1^{III}; \\ g\_{II}'(0) &= g\_{III}'(0): & -\frac{BL}{D\_{\mathcal{S}}} \mathcal{C}\_2^{II} + \frac{\mathcal{Q}}{BL} &= \frac{-BL - \sqrt{B^2 L^2 + 4D\_{\mathcal{S}}P}}{2D\_{\mathcal{S}}} \mathcal{C}\_1^{III}. \end{aligned}$$

Substituting the values of *CI I* 1 and *CI I* 2 into these equations allows obtaining the following implicit expression for *L*:

$$[1 - \frac{Q}{B} + \frac{QD\_{\mathcal{g}}}{B^2 L^2} [1 - e^{-BL^2/D\_{\mathcal{g}}}] [1 - \frac{2}{1 + \sqrt{1 + 4D\_{\mathcal{g}} P/(B^2 L^2)}}] = g\_{cr}. \tag{A1}$$

#### **References**

c 2020 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Optimization of Dose Fractionation for Radiotherapy of a Solid Tumor with Account of Oxygen Effect and Proliferative Heterogeneity**

**Maxim Kuznetsov 1,2,\* and Andrey Kolobov 1**

Received: 26 June 2020; Accepted: 18 July 2020; Published: 22 July 2020

**Abstract:** A spatially-distributed continuous mathematical model of solid tumor growth and treatment by fractionated radiotherapy is presented. The model explicitly accounts for three time and space-dependent factors that influence the efficiency of radiotherapy fractionation schemes—tumor cell repopulation, reoxygenation and redistribution of proliferative states. A special algorithm is developed, aimed at finding the fractionation schemes that provide increased tumor cure probability under the constraints of maximum normal tissue damage and maximum fractional dose. The optimization procedure is performed for varied radiosensitivity of tumor cells under the values of model parameters, corresponding to different degrees of tumor malignancy. The resulting optimized schemes consist of two stages. The first stages are aimed to increase the radiosensitivity of the tumor cells, remaining after their end, sparing the caused normal tissue damage. This allows to increase the doses during the second stages and thus take advantage of the obtained increased radiosensitivity. Such method leads to significant expansions in the curative ranges of the values of tumor radiosensitivity parameters. Overall, the results of this study represent the theoretical proof of concept that non-uniform radiotherapy fractionation schemes may be considerably more effective that uniform ones, due to the time and space-dependent effects.

**Keywords:** mathematical oncology; spatially-distributed modeling; reaction-diffusion-convection equations; computer experiment; gradient descent

**MSC:** 35K57; 35Q92; 92C05

#### **1. Introduction**

Approximately half of the patients, diagnosed with cancer, undergo radiotherapy (RT) [1]. The effect of irradiation on cancer and normal cells can be mathematically expressed via classical linear-quadratic model, which is known to fit experimental data well in a wide range of clinical parameters [2]. According to this model, the fraction of cells, which survive after a single radiation dose *D*, can be estimated as

$$S(D) = e^{-aD - \beta D^2} \, , \tag{1}$$

where *α* and *β* are radiosensitivity parameters of cells. Usually cancer cells have higher values of linear radiosensitivity parameter *α* that corresponding normal cells. However, normal tissues as a rule have greater *α*/*β* ratio, which restricts the use of high radiation doses [3]. One option to reduce normal tissue damage is to concentrate the radiation dose within the tumor mass. However, such option carries significant risks even for tumors with clear boundaries, due to the leakage radiation [4]. Other option to spare normal tissues is to fractionate the total dose, that is, to divide it into much smaller fractions, administered over a period of several weeks. The efficiency of a fractionation scheme depends on the effects that are widely referred to as the four "R"s of radiotherapy [5]. Two of the effects—reoxygenation and redistribution of cell cycle—indicate that the radiosensitivity of a cell depends on the surrounding concentration of oxygen and on the current stage of its cell cycle. In particular, hypoxic and non-proliferating cells are more radioresistant [6]. The third effect is the repopulation of tumor cells that takes place between the irradiations. The fourth effect is the repair of sublethal damage, which can be neglected, unless the time interval between irradiations is as short as several hours [7].

The most typical RT fractionation schemes consist of fractions of 1.8 to 2.0 Gy, delivered once a day on weekdays within the period of several weeks [3]. However, different fractionation protocols were shown to lead to improvement in tumor cure and patient survival for most of the patients for some of the tumor types [8]. Importantly, in clinical practice optimization of RT fractionation is significantly complicated by a number of factors, including the great variability of tumors, belonging to the same type, and related ethical problems, associated with the fact that alternative protocols may worsen outcome for some of the patients. Of note, in the vast majority of the tested schemes the irradiation doses are distributed equally between the fractions. At that, the varied parameters of the schemes are the number of the fractions, the interval between them and the fractional dose, which are related through the constraint on total normal tissue damage.

Given the practical difficulties, mathematical modeling can be a powerful tool, providing insights into the problem of optimization of RT fractionation. Different approaches exist in this field that have their pros and cons. The use of non-spatially distributed phenomenological models of tumor growth, expressed in ordinary differential equations, often allows to obtain the globally optimal solutions for the considered problems via analytical methods [9–11]. However, these methods become complicated and even unsolvable under introduction of complex non-linear terms, aimed to account for timeand space-dependent effects [12] or for tumor-specific features [13]. In such cases, various heuristic approaches are used, ranging from direct comparison of the schemes [14] to more complex techniques like simulated annealing [13], which can not guarantee the global optimality of the solution. However, such methods can yield significant results, one of which has already been verified in a preclinical study [13].

Crucially, the use of ordinary differential equations can allow for only phenomenological consideration of the reoxygenation and redistribution of the cell cycle—the effects that result are the spatiotemporal variability of tumor cells' radiosensitivity [12,14]. Their explicit consideration is possible in spatially-distributed models, which can be divided in two types—a continuous model, expressed in partial differential equations (PDEs), and discrete models that usually treat every tumor cell as a separate agent but use PDEs for the consideration of dynamics of nutrients and other substances. However, to the best of our knowledge, the existing works on RT fractionation optimization that use continuous spatially-distributed models, only account for homogeneous and constant radiosensivity of tumor cells. That leads to the conclusion of optimality of standard radiation fractionation schemes [15,16]. Reoxygenation and redistribution of cell cycle can be straightforwardly incorporated into agent-based models. However, the complexity of such models and numerical costs of their simulations lead to practical impossibility of solving optimization tasks via them, at least under current level of computer technology. As a rule, in the corresponding works several fractionation schemes are compared directly, moreover, the considered numbers of tumor cells are several orders of magnitude less that the relevant numbers for the human tumors [17–19]. These factors significantly limit the usefulness of such models.

In this work, we present a spatially-distributed continuous mathematical model of solid tumor growth and treatment by fractionated RT that explicitly accounts for tumor cell repopulation, reoxygenation and redistribution of proliferative states. With the use of a specially-developed algorithm, we find the optimized fractionation schemes for varied radiosensitivity of tumor cells

under the values of model parameters that correspond to different degrees of tumor malignancy. The resulting schemes lead to significant expansions in the curative ranges of the values of tumor radiosensitivity parameters.

#### **2. Model**

#### *2.1. Equations for Tumor Growth*

The mathematical model of tumor growth, considered herein, was based on our models, previously used for the investigation of various aspects of tumor growth and treatment [20–23]. There were five variables in this version of the model, which were the functions of space and time coordinates, *r* and *t*: the density of tumor cells *n*(*r*, *t*), the density of normal cells *h*(*r*, *t*), the fraction of necrotic tissue *m*(*r*, *t*), the concentration of glucose *g*(*r*, *t*) and the concentration of oxygen *ω*(*r*, *t*). The main simplification of this version of the model was the absence of an explicit variable for the capillaries, from which the nutrients flow into the tissue. Instead it was assumed that their local density was proportional to the local fraction of the normal cells in the tissue.

The following set of equations governed the dynamics of the model variables under the absence of radiotherapy:

$$\begin{array}{lcl} \text{tumor cells} & \stackrel{\text{followtion}}{\text{do}} & \stackrel{\text{death}}{\text{do}} & \text{migration correction} \\ \text{normal cells} & \stackrel{\text{dilat}}{\text{do}} & \stackrel{\text{followtion}}{\text{g}} \\ \text{normal cells} & \stackrel{\text{dilat}}{\text{do}} & \stackrel{\text{convolution}}{\text{cowection}} \\ \text{reccctic tissue} & \stackrel{\text{followtion}}{\text{do}} & \stackrel{\text{collash}}{\text{f}} & \text{cowection} \\ \text{reccctic tissue} & \stackrel{\text{followoff}}{\text{g}} & \stackrel{\text{followoff}}{\text{cow consumption}} & \text{cow consumption} \\ \text{glucose} & \stackrel{\text{followoff}}{\text{g}} & \stackrel{\text{conumption}}{\text{g}} & [\text{Q}^{\text{a}}\_{\text{f}} \text{in} + \text{Q}^{\text{g}}\_{\text{f}}\text{h}] \stackrel{\text{diffusion}}{\text{g}} & \stackrel{\text{diffusion}}{\text{cowection}} \\ \text{logover} & \stackrel{\text{followoff}}{\text{g}} & \stackrel{\text{followoff}}{\text{g}} & [\text{S}\_{\text{a}} \mu\_{\text{f}} + \text{Q}^{\text{g}}\_{\text{f}}] \stackrel{\text{diffusion}}{\text{g}} \\ \text{covogen} & \stackrel{\text{followoff}}{\text{cowignment}} & \stackrel{\text{conmission}}{\text{coworption}} \\ & & \text{coworption} & \\ \text{where} & -[\{Q^{\text{a}}\_{\text{g}} \frac{\mathcal{S}}{\mathcal{S} + \mathcal{S}^{\text{a}}} + Q^{\text{ai}}\_{\text{g}} \frac{\mathcal{S}^{\text{a}}}{\mathcal{S} + \mathcal{S}^{\text{a}}}\} + Q^{\text{a}$$

#### 2.1.1. Dynamics of Cells and Necrotic Tissue

The term of tumor cells proliferation implied that the rate of this process was proportional to the rate of glucose consumption by tumor cells. This assumption was made on the basis that glucose is indispensable nutrient for biosynthesis [24]. Glucose is also a crucial energetic nutrient for tumor cells [25]; however, they are known to obtain energy under its depletion via multiple ways in order to increase their survival [26–28]. Therefore, the limiting nutrient in the model for the cell survival was oxygen, tumor cells being more or at least equally resistant to its depletion that the normal cells, that is, ≤ 1. The function of cell death rate *Mh*(*ω*) was chosen to be smooth, tending to its maximum value under the full absence of oxygen, and equal to zero for the levels of oxygen which exceed the critical value *ω*∗. Normal cells also died in the presence of tumor cells, which was introduced in the model to coarsely reflect two processes: the inability of normal cells to remain viable in acidic tumor microenvironment [29] and the degradation of capillary network inside the tumor [30].

Tumor cells were able to migrate throughout the tissue, which was governed by a diffusion-like term. Directed migration of tumor cells was neglected [31]. The convective terms described the bulk motion of the tissue elements, the velocity field **I** being determined by the dynamics of tumor cells. The drainage of necrotic tissue was neglected. Due to the assumption of the constancy of the total density of tumor cells, necrotic tissue, and normal cells, which was normalized to unity, the following expression for the gradient of **I** was obtained:

$$
\nabla \mathbf{I} = B n \frac{\mathcal{S}}{\mathcal{g} + \mathcal{g}^\*} + D\_n \Delta n. \tag{3}
$$

#### 2.1.2. Dynamics of Nutrients

Model dynamics of both nutrients accounted for the same physiological processes—inflow of nutrients from the capillary network into the tissue, their consumption by tumor and normal cells and their diffusion within the tissue, the latter being much faster for oxygen. The form of the terms for the inflow of oxygen and glucose differed due to significant distinctions in the mechanisms of their blood and transvascular transport. The inflow of glucose is governed primarily by the process of passive diffusion through the pores in the walls of capillaries [32]. Therefore, the rate of the glucose inflow was taken to be proportional to the difference in glucose concentrations in blood and in tissue, and to the capillaries density, which, as it was mentioned above, was assumed to be in linear relationship with the density of normal cells. Glucose concentration in blood was considered to be constant and was normalized to unity.

Oxygen, as lipid-soluble substance with low molecular weight, passes directly through the capillary walls and flows into the tissue at much greater rate that glucose. Oxygen levels in arterial and venous blood differ more than twice even under normal conditions [33], which implies that its blood concentration should not be treated as constant. Moreover, the inflow of oxygen into the tissue is not proportional to the difference of its concentrations in capillary blood and tissue due to the complicated blood transport of oxygen, which molecules are carried in blood in two forms—bound to hemoglobin and unbound from it. Overall, the used term for the oxygen inflow assumed that the rate of this process is proportional to the difference between the fraction of oxygen-saturated hemoglobin under two values of unbound oxygen concentration—the one in arterial blood, which enters the capillaries, and the one in tissue. The function *S*(*ω*) represented oxygen-hemoglobin dissociation curve, the form of which is well-known in physiology [34]. For more detailed explanation of the assumptions, which underpinned the term of oxygen inflow, we refer the readers to our previous work [35].