Subdiffusions on circular branching structures

Commun Nonlinear Sci Numer Simulat 77 (2019) 225–238

Contents lists available at ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Research paper

Subdiffusions on circular branching structures Yu Fan a, Lin Liu a, Liancun Zheng a,∗, Xiaoliang Li b a b

School of Mathematics and Physics, University of Science and Technology Beijing, Beijing, 100083, China School of Mathematical Sciences, Beijing Normal university, 100875, China

a r t i c l e

i n f o

Article history: Received 14 December 2017 Revised 11 January 2019 Accepted 26 April 2019 Available online 3 May 2019 Keywords: Circular branching structure Subdiffusion

a b s t r a c t Anomalous diffusions on classical structures, like the comb structure where the branches are perpendicular to the straight backbone, have been studied by many researchers. However, in many practical problems, the anomalous diffusions are associated with a region surrounded by a boundary, and the research about such a kind of branching structures is still deficient. In this article, a novel model, namely, a circular backbone with radial branches, is presented to describe the anomalous diffusion on more complex regions, which has the potential of modeling many problems in medical science and nature, for example, angiogenic vasculatures induced by tumor. Three types of diffusion in circular branching structures are discussed, including two-way diffusion, inward diffusion, and outward diffusion. A significant advantage of these models is the solvability: The Fokker– Planck equations could be solved by variable separation, the Laplace transformation and its numerical inversion. The angular initial time mean squared displacement (MSD) is ob1 tained, and the results show that angular MSD is proportional to t 2 as t → 0, indicating that these types of diffusion are subdiffusion. Moreover, a numerical simulation of inward model is studied and the results are in good agreement with our analytical solution. Finally, radial MSD is calculated numerically, which shows that the radial diffusion is normal diffusion in short time, and the sink boundary conditions have a different behavior from the case with the reflective boundary conditions. © 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license. (http://creativecommons.org/licenses/by-nc-nd/4.0/)

1. Introduction There are many branching structures in nature [1], such as the spiny dendrites, the colonial development of bacteria, dendritic growth and so on. The diffusion phenomena in branching structures are of significant research as it can provide a deep understanding of nature. Many researchers have considered anomalous diffusions in classical comb structures (a well-known branching structure) and relevant applications. Arkhincheev and Baskin [2] considered the diffusion and drift of particles in percolation clusters 1

with comb structure, and MSD was obtained with < x2 (t ) >∼ t 2 . Fractal grid comb was introduced to describe the solvents diffusion in a thin porous film [3]. Baskin and Lomin [4] applied a comb model to tumor development and they proposed an explanation how malignant neoplasm cells can appear at an arbitrary distance from the primary tumor. Then a fractional, in time, tumor development was introduced into comb [5], and solid and diffusive cancers were discussed. Moreover, the tumor ∗

Corresponding author. E-mail addresses: [email protected], [email protected] (L. Zheng).

https://doi.org/10.1016/j.cnsns.2019.04.027 1007-5704/© 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license. (http://creativecommons.org/licenses/by-nc-nd/4.0/)

226

Y. Fan, L. Liu and L. Zheng et al. / Commun Nonlinear Sci Numer Simulat 77 (2019) 225–238

Fig. 1. A circular branching structure.

development through a lymphatic net as a kind of fractional transport of cells was suggested, and it showed that the main size of the growth is determined by subdiffusion, and that the metaphases is a result of superdiffusion [6]. The diffusion phenomena of calcium in dendrites [7,8] was investigated and an analysis of anomalous diffusion in spiny dendrites [9] was presented. Cattaneo–Christov flux and time were introduced to discuss the anomalous diffusion on comblike structures [10], so were spatial fractional derivatives [11]. However, many problems have more complex branching structures, and theoretical discussions of diffusion in such branching structures are still deficient. Here we take the angiogenic vasculature induced by tumor as example. When a vascular tumor is larger than 2 mm, the vasculature is no longer able to provide nutrients and oxygen to support its growth [12]. This will trigger the release of angiogenic factors, which will finally lead to the appearance of a finger-like capillary network around the tumor [13–20]. This process is called angiogenesis. The network is interconnected around the tumor and disjoint far from the tumor. If we consider this problem in 2D, the network interconnected around a circular tumor could be understood as a circumference (in our subsequent work it will be called a backbone), and the capillary relatively far away from the tumor can be seen as the branches. To treat the tumor, the injection of drugs is needed, and the concentration of drugs around the tumor would have significant influence on the treatment effect [13]. Chaplain et al. [21] investigated a discrete mathematical model of tumor-induced angiogenesis and the drug delivery to the tumor, namely, the diffusion in such a branching structure, by numerical simulation. Here, we use a theoretical model, which has a circular backbone with many secondary branches perpendicular to the backbone to describe diffusion in this kind of structures, which has never been reported before, to the best of our knowledge. Three cases of anomalous diffusions in circular branching structures, i.e., two-way diffusion, outward diffusion with infinite length branches and inward diffusion with branches of finite length l, are studied. We obtain the solution by performing the separation of variables, Laplace transformation and numerical inversion. Our results show that such a diffusion 1

phenomenon is a kind of subdiffusion along the circumference satisfying < θ 2 (t ) >∼ t 2 (t → 0 ), and normal diffusion along the branches satisfying < μ2 (t) > ∼ t, (t → 0). 2. Circular branching structures A special characteristic of diffusion in such circular branching structures is that the tangential current is absent as r = r0 , where r0 is the radial position where all particles were initially located. The diffusion on these structures is described by a 2D probability distribution function (PDF) P =P (t, r, θ ). We assume that all the particles are initially set at the position r = r0 and θ = 0, as shown in Fig. 1. P =P (t, r, θ ) is the probability of finding a particle at time t at the position r, θ in polar coordinates, implying that the particle has diffused to the backbone at a distance |θ r| and has diffused along the secondary branch at a distance |r − r0 |. 2.1. Two-way diffusion First, we consider the case of two-way diffusion where the particles could move in two directions: both inward and outward. Its currents are

Jr = −Dr

∂P r0 ∂ P , J = −Dθ δ ( r − r0 ), ∂r θ r ∂θ

Y. Fan, L. Liu and L. Zheng et al. / Commun Nonlinear Sci Numer Simulat 77 (2019) 225–238

227

where P =P (t, r, θ ) is the distribution function. For simplicity, we assume that Dr = Dθ = D. The Fokker–Planck equation, by an infinitesimal analysis, in polar coordinates, can be written as:

  ∂P 1 ∂P ∂ 2 P r0 ∂ 2 P =D + + δ ( r − r ) , 0 ∂t r ∂r ∂ r2 r2 ∂θ 2

(1)

where r ∈ (0, +∞ ). The initial condition

P (0, r, θ ) =

1 δ ( r − r0 )δ ( θ ) r0

(2)

means that all the particles are distributed at r = r0 and θ = 0 at the beginning, subject to the boundary conditions

0 ≤ P (t, 0, θ ) < ∞, P (t, ∞, θ ) = 0.

(3)

There is also a periodic condition

P (t, r, θ ) = P (t, r, θ + 2π ).

(4)

Because of the symmetry of the structure and the initial condition about x-axis, Jθ = 0 at θ = ±π , we have

∂P ∂P (t, r, π ) = (t, r, −π ) = 0. ∂θ ∂θ

(5)

Let t ∗ = Dt/r02 , P ∗ (t ∗ , r ∗ , θ ) = P (t, r, θ )r02 , and r ∗ = r/r0 . Then we have

∂ P∗ 1 ∂ P∗ ∂ 2 P∗ ∂ 2 P∗ ∗ = ∗ ∗ + + δ ( r − 1 ), ∗ ∗ 2 ∂t r ∂r ∂ (r ) ∂θ 2

(6)

P ∗ (0, r ∗ , θ ) = r02 P = δ (r ∗ − 1 )δ (θ ),

(7)

0 ≤ P ∗ (t ∗ , 0, θ ) < ∞, P ∗ (t ∗ , ∞, θ ) = 0,

(8)

P ∗ (t ∗ , r ∗ , θ ) = P ∗ (t ∗ , r ∗ , θ + 2π ),

(9)

and

∂ P∗ ∗ ∗ ∂ P∗ ∗ ∗ (t , r , π ) = (t , r , −π ) = 0. ∂θ ∂θ

(10)

For simplicity, we still denote P∗ as P, t∗ as t, and r∗ as r, respectively. Applying the Laplace transformation L[P (t, r, θ )] = P (s, r, θ ), Eqs. (6)–(10) become

sP −

∂ 2P 1 ∂ P ∂ 2P − − δ ( r − 1 ) = δ ( r − 1 )δ ( θ ), ∂ r2 r ∂ r ∂θ 2

(11)

0 ≤ P (s, 0, θ ) < ∞, P (s, +∞, θ ) = 0,

(12)

P (s, r, θ ) = P (s, r, θ + 2π ),

(13)

∂P ∂P (s, r, π ) = (s, r, −π ) = 0. ∂θ ∂θ

(14)

We assume that P (s, r, θ ) = g(s, θ ) f (s, r ) and substitute it into Eq. (11), yielding



g sf −

1∂f − r ∂r

∂2 f ∂ r2



− f (s, r )

∂ 2g δ ( r − 1 ) = δ ( r − 1 )δ ( θ ). ∂θ 2

(15)

Eq. (15) implies (see Appendix A)

sf −

1∂f − r ∂r

∂2 f = δ ( r − 1 ). ∂ r2

(16)

Therefore, we obtain the following two subproblems:



s f − 1r ∂∂ rf − ∂∂ r2f = δ (r − 1 ), f (s, ∞ ) = 0, | f (s, 0 )| < ∞, 2

(17)

228

Y. Fan, L. Liu and L. Zheng et al. / Commun Nonlinear Sci Numer Simulat 77 (2019) 225–238

Fig. 2. The probability concentration field of two-way model at t=1 (r ≤ 2).

and



∂ g g − f (s, 1 ) ∂θ 2 = δ ( θ ), 2

(18)

∂g ∂g ∂θ (s, π ) = ∂θ (s, −π ) = 0,

which are akin to Green’s function problems, and their solutions are obtained as

f = C1 I0

√ 

sr (1 − η (r − 1 )) + η (r − 1 )C2 K0

√ 

sr ,

g(s, θ ) = [a1 exp(λ|θ | ) + a2 exp(−λ|θ | )], where

(19)

√ K0 ( s ) 1 C1 = √ √ √ √ √ , K0 ( s )I1 ( s ) + I0 ( s )K1 ( s ) s √ I0 ( s ) 1 C2 = √ √ √ √ √ , K0 ( s )I1 ( s ) + I0 ( s )K1 ( s ) s

 λ= a1 =

1 , f (s, 1 )

λ

2(exp(2λπ ) − 1 )

, a2 =

λ exp(2λπ ) , 2(exp(2λπ ) − 1 )

and In is n-th order modified Bessel function of first kind, while Kn is of the second kind, and η (r − 1 ) is a heaviside function. Combining these two solutions, we find the solution in the form P (s, r, θ ) = g(s, θ ) f (s, r ). The graphics of concentration field is depicted in Fig. 2. Furthermore, its angular MSD is obtained as follows,

 

π 2π exp (λπ ) θ 2 P (t, r, θ )dθ L−1 f (s, r ) 0 θ 2 g(s, θ )dθ L−1 f (s, r ) λ12 − λ(exp (2λπ )−1)

= < θ (r, t ) >= = . π

π 2 0 P (t, r, θ )dθ L−1 f (s, r ) 0 g(s, θ )dθ L−1 f (2s,r ) 2

2

π 0

(20)

Y. Fan, L. Liu and L. Zheng et al. / Commun Nonlinear Sci Numer Simulat 77 (2019) 225–238

229

Fig. 3. The time dependence of angular MSD on the circular structure in a two-way diffusion model (in long time and r = 1).

Fig. 4. The growth pattern of angular MSD on t1/2 (in short time).

Especially, on the backbone, the angular MSD is



< θ ( 1, t ) >=< θ 2

2

π exp (λπ ) L−1 λ12 λ12 − λ(2exp 1 (2λπ )−1) (t ) >= L−1 2λ2



.

(21)

Compared to the angular MSD on the circumference r = r0 = 1, the angular MSD on the circumference r = 1 is more important, because the particles could only move angularly along the backbone. The angular MSD < θ 2 (t) > on the structure can be calculated by Laplace numerical inversion [22], and the result is plotted in Fig. 3. Through the observation from Fig. 3, we find that the angular MSD has a linear growth pattern in short time, which is shown in the following picture in a smaller time scale and with different r. Moreover, Fig. 4 shows that the growth patterns of angular MSD on different circumferences are the same if their circumferences have the same distance to r = 1.

230

Y. Fan, L. Liu and L. Zheng et al. / Commun Nonlinear Sci Numer Simulat 77 (2019) 225–238

2.2. Outward diffusion Now, we consider an outward diffusion whose spatial domain is r ≥ 1. We face a challenge that the particles are distributed on the boundary at initial time, which is hard to describe by initial condition on open interval. Therefore, we extend the spatial domain from r ≥ 1 to r ≥ a, where 0 < a < 1, and add a source ∂ P/∂ t |(t,1+ ,θ ) (1 − η (r − 1 )) on the extended region (inspired by the electric image method). Then, the diffusion equation is turned into

  ∂P 1 ∂ ∂P ∂ 2P ∂P − r − δ (r − 1 ) = (t, 1+ , θ )(1 − η (r − 1 )), ∂t r ∂ r ∂ r ∂t ∂θ 2

(22)

with the initial distribution P (0, r, θ ) = δ (r − 1 )δ (θ ), and boundary conditions

P (t, +∞, θ ) = 0,

∂P (t, a, θ ) = 0. ∂r

(23)

If we extend the spatial domain in such a way, there would be a particle current flowing from circumference to the interior. We have to add a source (the right hand of Eq. (22)) which could eliminate the concentration gradient inside the circumference and prevent the particles from flowing inward. Performing the Laplace transform on Eq. (22), we obtain



sP −

1 ∂ ∂P r r ∂r ∂r



−

∂ 2P δ (r − 1 ) = sP (s, 1+ , θ )(1 − η (r − 1 )) + δ (r − 1 )δ (θ ). ∂θ 2

Taking P (s, r, θ ) = f (s, r )g(s, θ ), we obtain





1 ∂ ∂f g sf − r r ∂r ∂r Let



sf −

1 ∂ ∂f r r ∂r ∂r





− s f (s, 1 )(1 − η (r − 1 )) − +

∂ 2g f (s, 1 )δ (r − 1 ) = δ (r − 1 )δ (θ ). ∂θ 2

 − s f (s, 1+ )(1 − η (r − 1 )) = δ (r − 1 ).

(24)

We consider this equation in two parts. Obviously, for r > 1, the solution is f = CK0 and then we get the derivative of f as

√ 

sr . For r < 1, we have f = f (s, 1+ ),

∂f 0, a√< r < √1, = − C sK1 ( sr ), r > 1. ∂r

Multiplying the both sides of Eq. (24) by r, we obtain

  ∂ ∂f sr f − r − sr f (s, 1+ )(1 − η (r − 1 )) = δ (r − 1 ). ∂r ∂r

Integrating Eq. (25) from 1− to 1+ with respect to r, we obtain − ∂∂ rf |1+ = 1 and C = 1−

⎧ √ K ( sr ) ⎪ ⎨ 0 √ √1s , r ≥ 1, K1 (√ s ) f = ⎪ ⎩ K0 (√s ) √1 , a ≤ r < 1. K1 ( s ) s

(25) 1√ √ 1 . K1 ( s ) s

Therefore, we have

Utilizing the solution of f (r ≥ 1), we obtain the solution g which has the same form as Eq. (19). The concentration field of an outward model is depicted in Fig. 2. Furthermore, we can calculate the MSD in a similar way as in the two-way diffusion case. In fact, when calculating, we mainly focus on r ∈ (1, +∞ ), as the inside domain is only auxiliary.(Fig. 5) 2.3. Inward diffusion with branches of finite length l Similarly, if the branches only exist inside the circumference with a length l, we need to extend the spatial domain from r ≤ 1 to r ≤ a (a > 1). Correspondingly, the diffusion equation is

  ∂P 1 ∂ ∂P ∂ 2P ∂P − r − δ (r − 1 ) = (t, 1− , θ )η (r − 1 ), ∂t r ∂ r ∂ r ∂t ∂θ 2

with the same initial conditions as before. But the boundary condition at the end of the branch is turned into ∂∂Pr |r=1−l = 0 because the length of branches is l. The other boundary condition is ∂∂Pr |r=a = 0. We obtain the solution of f in a similar way, as

f =

√ √ B1 I0 ( sr ) + B2 K0 ( sr ), 1 − l ≤ r < 1, √ √ B1 I0 ( s ) + B2 K0 ( s ), 1 < r < a.

Y. Fan, L. Liu and L. Zheng et al. / Commun Nonlinear Sci Numer Simulat 77 (2019) 225–238

231

Fig. 5. The probability concentration field of the outward model at t=1 (r ≤ 2).

where

√ √ K1 ( s(1 − l ))/ s √ √ √ , K1 ( s(1 − l ))I1 ( s ) − I1 ( s(1 − l ))K1 ( s ) √ √ I1 ( s(1 − l ))/ s B2 = √ √ √ √ . K1 ( s(1 − l ))I1 ( s ) − I1 ( s(1 − l ))K1 ( s ) B1 =

√

Similarly here, the solution of g has the same form as Eq. (19). Utilizing the solution of f and g, we can calculate the MSD in a similar way as in the two-way diffusion case. Here, the outside domain is auxiliary. The concentration field of the inward model (l = 1) is depicted in Fig. 6. Let θ = 0. We calculate the probability distribution with respect to r. Noting that the total probability is  P (t, r, θ )rdrdθ = 1, and the area of the circle ring is S = π (2l − l 2 ), as a result, we find that the probability at each point converges to 1/π (2l − l 2 ) as t → ∞. Take l → 1 for example. The circular ring becomes a disk with an area π , and the probability at each point converges to 1/π = 0.3183, which coincides with our calculation result well and is shown in Fig. 7. 2.4. Growth pattern of MSD near initial time To discuss the growth pattern of MSD, we need to introduce two important propositions. 1

Proposition 1. If L[ f (t )] = f (s ), and f (s ) ∼ s− 2 ,

df (s ) ds

3

1

∼ s− 2 as s → ∞, then f (t ) ∼ t 2 as t → 0.

Proof. Denote g(t ) = t f (t ) and L[g(t )] = g(s ) = − dfds(s ) ∼ s− 2 . Then we have limt→0 g(t ) = lims→∞ sg(s ) = 0. Considering the fractional derivative and integral of g: 3



−1



1

L Dt 2 g(t ) = s− 2 g(s ) ∼ s−2 .

(26)

According to initial value theorem, we have −1

1

lim Dt 2 g(t ) = lim ss− 2 g(s ) = 0, s→∞

t→0

and



1



1



−1

L Dt2 g(t ) = s 2 g(s ) − Dt 2 g(t )

 1 |t=0 = s 2 g(s ) ∼ s−1 .

(27)

(28)

232

Y. Fan, L. Liu and L. Zheng et al. / Commun Nonlinear Sci Numer Simulat 77 (2019) 225–238

Fig. 6. The probability concentration field of an inward model at t=1 (r ≤ 2).

Fig. 7. Probability distribution with respect to θ on the circumference r = 1 in the inward diffusion model (l = 1).

So, estimating initial value from the Laplace transform properties again, we have 1

1

lim Dt2 g(t ) = lim ss 2 g(s ) = c = 0. ss→∞

t→0

(29)

For all ε > 0, there exists a δ > 0, when 0 < t < δ , such that 1

c − ε < Dt2 g(t ) < c + ε .

(30)

Y. Fan, L. Liu and L. Zheng et al. / Commun Nonlinear Sci Numer Simulat 77 (2019) 225–238

233

Considering the 1/2 order fractional integral of above equation, −1

−1

1

−1

(c − ε )Dt 2 1 < Dt 2 Dt2 g(t ) < (c + ε )Dt 2 1, namely

√ 2 t

(31)

√ 2 t

(c − ε ) √ < g(t ) < (c + ε ) √ , π π we have

−ε <

√

π f (t ) 1

2t − 2

(32)

− c < ε,

(33)

which yields that 1

f (t ) ∼ t − 2 , (t → 0 ).

(34) 

Proposition 2. If L[ f (t )] = f (s ) and f (s ) ∼ s−1 as s → ∞, we have f(t) → const = 0 as t → 0. Proof. From the Laplace transform properties, we have

lim f (t ) = lim s f (s ) = const

(35)

s→∞

t→0

 With these propositions, we can calculate the order of numerator and denominator of angular MSD at an initial time by letting g(t ) = t α . Furthermore, we can obtain the growth pattern of angular MSD near the initial time. For three cases above, due to asymptotic form of bessel functions:

Im (x ) ∼

ex π e−x √ , Km (x ) ∼ √ , (x → ∞ ), 2 x 2 x

we have

1 1 ∼ s− 2 , 2λ2 1



1

λ2 λ2

(36)

2π exp (λπ ) − λ(exp (2λπ ) − 1 )

Then, we have

 π 0

P (t, 1, θ )dθ ∼ t − 2 , 1

 π 0

 ∼ s−1 (s → ∞ ).

(37)

θ 2 12 P (t, 1, θ )dθ ∼ 1 (t → 0 ). 1

Therefore, we obtain the growth pattern of angular MSD at initial time in three cases < x2 (t ) >∼ t 2 (t → 0). 2.5. Proof of the Existence of Asymptote of angular MSD in the inward diffusion model Recall the definition of angular MSD Eq. (20). Note that probability P(t, 1, θ ) ≥ 0, P(t, 1, θ ) and θ 2 are continuous with respect to θ on the interval [0, π ]. Based on the mean value theorem of integrals, we have t→+∞

π

θ 2 12 P (t, 1, θ )dθ t→+∞ 0 P (t, 1, θ )d θ π P (t, 1, θ )dθ = θ02 lim 0π = θ02 , t→+∞ P ( t, 1 , θ ) d θ 0

lim < θ 2 (t ) > = lim

0

π

where θ 0 ∈ (0, π ). So the asymptotic value of angular mean squared displacement on the backbone exists. Combining with our steady solution of inward model, we know that in the inward model θ0 = √π , because P (t, 1, θ ) = π1 as t → +∞. 3

2.6. MSD along the radial direction The angular MSD has been discussed above. In this section, we will further discuss the MSD along the radial direction, which indicates the characteristic of the diffusion on the teeth. Here, we take the outward diffusion (along the teeth with finite length l), for example, to illustrate our definition.

234

Y. Fan, L. Liu and L. Zheng et al. / Commun Nonlinear Sci Numer Simulat 77 (2019) 225–238

Fig. 8. The local coordinate system at (r, θ ) = (1, θ ).

As Fig. 8 shows, we introduce a local coordinate system, whose x -axis is parallel to the radial direction, y -axis is perpendicular to the x -axis and origin coincides with (1, θ ). The concentration along the teeth θ = θ0 is Pθ˜0 (t, x , 0 ) = P (t, 1 + x , θ ), and the relation between x and r is x = r − 1. Then the radial MSD becomes

 y /2  b−1 

P˜θ0 (t, x , 0 )x2 dx dy

/2 0 μ(t, θ ) =  −y

y /2  b−1

 ≈

−y /2 0

i

P˜θ0 (t, x , 0 )dx dy

P (t, ri , θ )(ri − 1 ) r  = i P (t, ri , θ ) r 2



i

 b−1

=

P˜θ0 (t, x , 0 )x2 dx =  b−1 P˜θ (t, x , 0 )dx

0

0

0

b 1

P (t, r, θ )(r − 1 )2 dr b 1 P (t, r, θ )dr

P (t, ri , θ )(ri − 1 )2  . i P (t, ri , θ )

(38)

This formula can be calculated numerically; it also holds for the inward diffusion model with teeth with a finite length 0 < l < 1. 2.7. Discussion on radial MSD Based on the definition of radial MSD, we consider the inward diffusion with teeth with length l and calculate it numerically (which can be derived similarly), and the results are shown in Figs. 9 and 10. Fig. 9 shows that in short time, the growth of radial MSD is linear, indicating a kind of normal diffusion. Moreover, the longer teeth means that the ends are farther away from the backbone, which is beneficial to the diffusion along the teeth. Fig. 10 depicts the relation between the radial MSD and the angle: The radial MSD on the branch of θ = 0 grows the fastest and on the branch of θ = π grows the slowest. Figs. 11and 12 show that, all the cases have the same characteristic that the radial MSD has a linear growth pattern in short time. But compared to the reflective boundary conditions, the farther of the end (sink boundary) from the backbone, the faster the particles diffuse. 3. Numerical simulation We now provide a numerical simulation of the inward diffusion model in order to verify our results and illustrate the physical point further (see Fig. 7), and the details are shown as follows. In our simulation, the finite volume method is adopted, and the nodes are set at the center of the subdomains. We employ Mr + 1 nodes along the radial direction and Mt nodes along the tangential direction, and set Mr = Mt = 50, t = 10−5 (the distance of time step). The corresponding discrete scheme is

     

t

r  k

r  k k k = + ri + Pi+1, j − Pi, j − ri − Pi, j − Pi−1, j , (2 ≤ i ≤ Mr ), 2 2 r i r 2  4 t  k Pi,k+1 = Pi,k j + P − P1k, j , (i = 1 ), j

r 2 2 , j Pi,k+1 j

Pi,k j

and for i = Mr + 1

PMk+1 = PMk r +1, j + r +1, j



2 t 1 − 2r

r −

r 2 4



 2 t  k  1  k 1 k k PMr , j − PMk r +1, j + PMr+1, j+1 − 2PMr+1 2 , j + PMr+1, j−1 . 2

r

r

θ r − 4

Y. Fan, L. Liu and L. Zheng et al. / Commun Nonlinear Sci Numer Simulat 77 (2019) 225–238

235

Fig. 9. The dependence of radial MSD in the inward diffusion model along teeth with length l on the length (θ = 0 ).

Fig. 10. The dependence of radial MSD in inward diffusion along teeth with length l = 0.9 on the angle.

The initial condition is chosen to be

PM1 r +1,1 =

2



θ r −

r 2

.

4

k Note that P is periodic in the angular direction, so Pi,M

t+1

k . → Pi,k1 , and Pi,k0 → Pi,M t

Moreover, the reference time of subdiffusion at the initial time can be obtained by our numerical simulation. Because of the effect of the boundary (along the angular direction), the characteristics of subdiffusion fades away. In our simulation, during one step, the particles can only diffuse from one unit on the circumference to another neighbouring unit. So for particles, it needs θπ/ t steps to arrive at the boundaries, and the reference time can be represented as

tre f = lim

θ → 0

t→0

π π r0 π r0 = lim = ,

θ / t θ →0 θ r0 / t v0

t→0

(39)

236

Y. Fan, L. Liu and L. Zheng et al. / Commun Nonlinear Sci Numer Simulat 77 (2019) 225–238

Fig. 11. The dependence of radial MSD in the inward diffusion model with sink boundary conditions on t (at r = 1 − l).

Fig. 12. The dependence of radial MSD in the outward diffusion model with sink boundary conditions on t (at r = 1 + l).

where v0 is the mean velocity of the particles on the circumference r = r0 .

4. Conclusion In summary, we have mainly investigated three cases of diffusion in a circular branching-like structure. The analytical solutions after Laplace transform are obtained, and concentration fields and then MSD are calculated numerically. What is more, a numerical simulation of an inward model is provided and compared with our results. In three cases, all MSDs are rapidly increasing with respect to time at first. Because of the closed nature of circumference, the growth of MSDs are limited and finally convergent to a constant (see Section 2.5). Therefore, it is impossible to find a growth pattern tα to describe MSD in the whole time.

Y. Fan, L. Liu and L. Zheng et al. / Commun Nonlinear Sci Numer Simulat 77 (2019) 225–238

237

However, we can consider the growth patterns of angular MSD at initial time when the diffusion has not been affected 1

by the boundary remarkably. The angular MSDs at initial time have the form of t 2 (see Section 2.4), which indicates that the anomalous diffusion on the circumference is subdiffusion. This characteristic coincides with the results obtained in the classical comb model in the literature [14]. These results also imply that the presence of secondary radial branches plays a similar role in the diffusion as ”dead ends” in the reference, which leads to the anomalous diffusion whose MSDs are 1 proportional to t 2 as t → 0. In addition to the angular MSD, we have also introduced a local coordinates to study the radial MSD. We find that the radial MSD grows linearly on the branch, which is a kind of normal diffusion, and the sink boundary conditions have an different effects to the reflective boundary conditions: the sink boundary conditions accelerate the diffusion along the branch while the reflective boundary conditions restrain the diffusion. Moreover, according to the probability of the concentration field of the three cases, we conclude that the existence of branches enhances the diffusion along the radial direction. The outward model can explain why the rate of drug delivery to the tumor is very low. On the contrary, the inward model can simulate the diffusion of drugs without outside capillary and the two-way model can mimic the diffusion phenomenon with both inside and outside capillary. Acknowledgments The work is supported by the National Natural Science Foundations of China (no. 11772046, 81870345), and China Postdoctoral Science Foundation (no. 2018T110042). During revising this paper, we obtained many enthusiastic help from Professor Goong Chen (Texas A&M University, College Station, Texas), who checked whole manuscript and offered many valuable suggestions. we wish to express our sincerest thanks to Professor Goong Chen’s assistance. Appendix A. Derivation of Eq. (16) 2 2 Due to s f − 1r ∂∂ rf − ∂∂ r2f = 0 as r = 1, we suppose that s f − 1r ∂∂ rf − ∂∂ r2f = k(s )δ (r − 1 ) (obviously k(s) = 0). If we assume that the solution of Eq. (16) is f and the solution of Eq. (18) is g. Then the solution of above equation is k(s)f. Correspondingly, Eq. (18) becomes



∂ g k(s )g − k(s ) f (s, 1 ) ∂θ 2 = δ ( θ ), 2

∂g ∂g ∂θ (s, π ) = ∂θ (s, −π ) = 0,

and its solution is g/k(s), then we have P = (k(s ) f )(g/k(s )) = f g, which implies that let k(s ) = 1 is reasonable. Appendix B. The discontinuity at the origin It is worth noting that, in the inward case when l → 1 or in the two-way case, if we let r → 0+ , we could find that limr→0+ P (s, r, θ ) = g(s, θ ) limr→0+ f (s, r ). For limr→0+ f (s, r ), however, it is not convergent to 0, which indicates that the distribution function is discontinuous at the origin. Now, we consider the angular current of a normal inward diffusion in ∂ P . When r tends to 0+ , the tangential current tends to infinity. If there is a non-zero polar coordinates first, Jθ = −Dθ 1r ∂θ concentration gradient in the tangential direction, the infinite tangential current will eliminate the gradient immediately. In our two-way model and inward diffusion model, however, the angular current is limited so strictly that the tangential current is absent around the origin, which results in the discontinuity at the origin. In fact, our models describe a circular branching diffusion on a plane without the origin, or from physical standpoint, the branches are relatively independent, and so they are disjoint at origin. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

Fleury V, Gouyet J-F, Léonetti M. Branching in nature. Springer-Verlag Berlin Heidelberg 2001; 2001. Arkhincheev VE, Baskin EM. Anomalous diffusion and drift in a comb model of percolation clusters, 300; 1991. p. 161–5. Shamiryan D, Baklanov MR, Lyons P, Beckx S, Boullart W, Maex K. Diffusion of solvents in thin porous films. Colloids Surf A 20 07;30 0:111–16. Baskin E, Iomin A. Superdiffusion on a comb structure. Phys Rev Lett 2004;93. Iomin A. Toy model of fractional transport of cancer cells due to self-entrapping. Phys Rev E 2006;73:1–5. Iomin A, Dorfman S, Dorfman L. On tumor development: fractional transport approach; 2004. Santamaria F, Wils S, Schutter ED, Augustine GJ. Anomalous diffusion in purkinje cell dendrites caused by spines. Neuron 2006;52:635–48. Korkotian E, Segal M. Spatially confined diffusion of calcium in dendrites of hippocampal neurons revealed by flash photolysis of caged calcium. Cell Calcium 2006;40:441–9. Iomin A, Méndez V. Reaction-subdiffusion front propagation in a comblike model of spiny dendrites. Phys Rev E 2013;88:1–7. Liu L. Fractional anomalous diffusion with Cattaneo–Christov flux effects in a comb-like structure. Appl Math Model 2016;40:6663–75. Liu L, Zheng L, Liu F, Zhang X. Anomalous convection diffusion and wave coupling transport of cells on comb frame with fractional Cattaneo–Christov flux. Commun Nonlinear Sci Numer Simul 2016;38:45–58. Folkman J. Tumor angiogenesis: therapeutic implications. N Engl J Med 1971;285:1182. Sanga S, Sinek J, Frieboes H. Mathematical modeling of cancer progression and response to chemotherapy. Expert Rev Anticancer Ther 2006;6:1361–76. Baskin E, Iomin A. Superdiffusion on a comb structure. Phys Rev Lett 2004;93:1–4. Lee DS, Rieger H, Bartha K. Flow correlated percolation during vascular remodeling in growing tumors. Phys Rev Lett 2006;96:1–4.

238

Y. Fan, L. Liu and L. Zheng et al. / Commun Nonlinear Sci Numer Simulat 77 (2019) 225–238

[16] Carmeliet P, Jain RK. Angiogenesis in cancer and other diseases. Nature 20 0 0;407:249–57. [17] Ausprunk DH, Folkman J. Migration and proliferation of endothelial cells in preformed and newly formed blood vessels during tumor angiogenesis. Microvasc Res 1977;14:53–65. [18] Sholley MM, Ferguson GP, Seibel HR, Montour JL, Wilson JD. Mechanisms of neovascularization. vascular sprouting can occur without proliferation of endothelial cells.. Lab Investig J Tech Methods Pathol 1984;51:624–34. [19] Gimbrone MA, Cotran RS, Leapman SB, Folkman J. Tumor growth and neovascularization: an experimental model using the rabbit cornea. J Natl Cancer Inst 1974;52:413–27. [20] Muthukkaruppan VR, Kubai L, Auerbach R. Tumor-induced neovascularization in the mouse eye. J Natl Cancer Inst 1982;69:699–708. [21] Chaplain MA, McDougall SR, Anderson AR. Mathematical modeling of tumor-induced angiogenesis. J Math Biol 2004;49:111–87. [22] Valsa J, Lubomír B. Approximate formulae for numerical inversion of Laplace transforms, 11; 1998. p. 153–66.