Fluctuation-induced dissipation in evolutionary dynamics

Tsung-Cheng Lu1, Yi-Ko Chen1, Hsiu-Hau Lin1 and Chun-Chung-Chen2

1Department of Physics, National Tsing Hua University, 30013 Hsinchu, Taiwan2 Institute of Physics, Academia Sinica, Nankang, Taipei 11529, Taiwan

Correspondence and requests for materials should be addressed toH.-H. L. (hsiuhau.lin@gmail.com) or C.-C. C (cjj@phys.sinica.edu.tw)

(Dated: October 4, 2014)

Biodiversity and extinction are central issues in evolution. Dynamical balance among differentspecies in ecosystems is often described by deterministic replicator equations with moderate success.However, fluctuations are inevitable, either caused by external environment or intrinsic random com-petitions in finite populations, and the evolutionary dynamics is stochastic in nature. Here we showthat, after appropriate coarse-graining, random fluctuations generate dissipation towards extinc-tion because the evolution trajectories in the phase space of all competing species possess positivecurvature. As a demonstrating example, we compare the fluctuation-induced dissipative dynamicsin Lotka-Volterra model with numerical simulations and find impressive agreement. Our finding isclosely related to the fluctuation-dissipation theorem in statistical mechanics but the marked differ-ence is the non-equilibrium essence of the generic evolutionary dynamics. As the evolving ecosystemsare far from equilibrium, the relation between fluctuations and dissipations is often complicated anddependent on microscopic details. It is thus remarkable that the generic positivity of the trajectorycurvature warrants dissipation arisen from the seemingly harmless fluctuations. The unexpecteddissipative dynamics is beyond the reach of conventional replicator equations and plays a crucialrole in investigating the biodiversity in ecosystems.

Biodiversity is commonly used to indicate the sta-bility of an ecosystem[1–3]. One of the central is-sues is to effectively promote the biodiversity while at-tracting more scientists’ attention from various fields[4–8]. The causes that threaten the biodiversity, for in-stance, climate change[4, 5], over-harvesting, habitatdestruction[6], and population mobility[7], are well stud-ied. Above those factors, Darwin’s theory of natural se-lection plays a crucial role in catalysis[9–12]. People arewarned to reduce these effects in order to maintain andreserve the nature’s biodiversity. Nevertheless, a naivereversed statement should be check, namely, would theecosystem be perfectly stable without any of these haz-ardous factors?

To put the discussions on firm ground, we can startwith the non-transitive rock-paper-scissors game[13–18],known as a paradigm to illustrate the species diver-sity. When three subpopulations interact in this non-transitive way, we expect that each species can in-vade another when its population is rare but becomesvulnerable to the other species when over populated.The non-hierarchical competition[19–22] gives rise to theendlessly spinning wheel of species chasing species andthe biodiversity of the ecosystem reaches a stable dy-namical balance. This cyclic evolutionary dynamicshas been found in plenty of ecosystems such as coralreef invertebrates[23], lizards in the inner Coast Rangeof California[24] and three strains of colicinogenic Es-cherichia coli [14, 25] in Petri dish. Although the os-cillatory solutions for the replicator equations capturethe main features, inclusion of mobility[7] or/and finite-population effects[21, 22] in the numerical simulationsalways jeopardizes the stable equilibrium and highlight

the importance of stochasticity in the evolution.Populations in an ecosystem are discrete integers.

Approximating these discrete populations by continu-ous variables inevitably introduces intrinsic fluctuations,which turn the evolutionary dynamics stochastic in na-ture. By extensive numerical simulations, we record howthe biodiversity of the ecosystem dissipates from theseintrinsic fluctuations. In addition to the irregular devi-ations at the short-time scale, slowly but surely, dissi-pation induced from intrinsic fluctuations is derived andlead to an extinction time proportional to the populationsize[26, 27]. Our findings can be elegantly summarizedin three steps: discreteness induces fluctuations, fluctua-tions spawn dissipations and dissipative dynamics leadsto extinction.

RESULTS

Intrinsic fluctuations. We start with the simplest casewith two species A and B, whose population numbers aredenoted by x1 and x2 respectively. Because the change ofthe population number of species is always discrete, thephase space spanned by x1 and x2 has internal uniformgrid structure, where the system can only evolve on gridpoints. For a given state of system on a grid point in asmall time step, the system can either hop to adjacentpoint or stay at the same point. Therefore, the systemis stochastic in nature and has to be described by proba-bility which is determined by competing relation amongspecies. As long as the total population N is large, withan appropriate coarse graining, we can derive effectivereplicator equations with details found in the Methods

arX

iv:1

411.

6473

v1 [

q-bi

o.PE

] 2

4 N

ov 2

014

2

X1

X2

Phase space structure

FIG. 1: Evolutionary contours in the phase space and theenlarged segment with local coordinates.

section,

dxidt

= fixi + ξi i = (1, 2) (1)

where fi stands for the fitness which governs the deter-ministic dynamics in the system, and ξi is the intrinsicnoise arisen from the discreteness of population. Notethat ξi approximates to Gaussian noise under the large-N limit. If deterministic dynamics of system driven byfixi predicts a closed trajectory in the phase space, thereexists a constant of motion which is determined by theinitial condition (x1(0), x2(0)), and thus the system canevolve on different contours depending on different ini-tial condition(Fig.1). However, system will deviate fromthe original orbit due the random noises. To capturethe influence caused by the noises, we can define stabil-ity indicator χ, which is a conserved quantity within thedeterministic dynamics, and let χ be the maximum atthe fixed point predicted by the deterministic dynamics.With this indicator, we can measure the effect of fluctu-ations by studying the time evolution of χ.

To realize the influence of noise, let us focus on the sys-tem at point A on the orbit as shown in Fig. 1. Duringa small time interval ∆t, the system will have a stochas-tic change in the direction of x1 and x2 caused by ξ1and ξ2 respectively. If we mark all the probable pointswith the same probability that the system can arrive in∆t, then we will get a closed orbit Γ which shows re-flectional symmetry with the vertical and horizontal axiswith respect to point A due to the isotropy of noise underlarge-N limit. Due to the curvature of equi-χ contour, thearea surrounded by blue line corresponding to χ > χ0 isgreater than the area surrounded by red line correspond-ing to χ < χ0. Therefore, the expectation value of thechange of χ will be negative(i.e.〈∆χ〉 < 0), which indi-cates that the dissipation of χ depends on the curvatureof equi-χ contours. Because the ratio of these two ar-eas will increase when the curvature increases, we canimagine that the magnitude of dissipation of χ and lo-cal curvature is positive correlated. Besides, for a givendistance that system can travel, ∆χ is larger in the areawith denser equi-χ contours, which is a reflection of themagnitude of gradient |∇χ|. Under large-N limit, for the

FIG. 2: (a)∆χ calculated by algorism in phase space witheach species number ranges from 20 to 90. (b)∆χ calculatedby the theory in phase space with each species number rangesfrom 20 to 90.

case that the standard deviations of Gaussian noises σin both vertical and horizontal directions are the same,the probable hopping points correspond to equal proba-bility becomes a circle centered at point A. If we take thesmall curvature approximation κ � 1

σ and assume thatlocal curvature of contours κ and |∇χ| are constant, theformula for 〈∆χ〉 can be derived explicitly:

〈∆χ〉 = −1

2κ|∇χ|σ2 (2)

The most important feature in Eq.(2) is that the mag-nitude of dissipation of χ has both linear dependanceon curvature and gradient, which reveals that the dissi-pation is determined by the geometric structure of con-tours in the phase space. Additionally, the magnitudeis proportional to σ2, which reflects the fact the dissipa-tion occurs due to the fluctuations in the system. Fromthe argument above, for those system with convex con-tours in phase space, 〈∆χ〉 is always negative, and thusthe system is inclined to move to outer and outer con-tours. When the system arrive at the boundaries, onlyone surviving species is left, revealing the destruction ofbiodiversity. Practically speaking, because any kind ofspecies cannot survive without the interaction with otherspecies, the ecosystem is destined to go to extinction.

The above method can be easily generalized to N-species stochastic system, where N can be any finite pos-itive integer. In N-dimensional phase space, if the deter-ministic dynamics predicts a closed convex orbit, where

3

we can define a biodiversity indicator χ for a orbit, thestructure of phase space can be characterized by a hy-persurface for different value of χ. During a small timeinterval ∆t, for a given state on a contour, we can plot ahypersurface Γ which is the set of equally probable pointsthat the system can hop into. Due to the convexity ofequi-χ hypersurface, the volume closed by the hypersur-face Γ is divided to two parts with different volumes. Inconsequence, the system is inclined to go to outer hyper-surface, implying the destined extinction of the system.Therefore, for any dimensional phase space where deter-ministic driving force predicts a closed topological struc-ture, the seemly harmless fluctuation can always inducethe dissipation, and thus drive the system to extinction.Lotka-Volterra Model. We present Lotka-Volterramodel as a demonstration of our method. The Lotka-Volterra model expresses the dynamics of the prey andthe predator in an ecosystem. Eqs.(7) are the replicatorequation, with with x and y denoting species populationof X(predator) and Y(prey) respectively. Besides the in-teraction part, there are natural birth of X and naturaldeath of Y using coefficient a and c in the mechanism.

dx

dt= −ax+ bxy

dy

dt= cy − dxy (3)

It is clear that the populations of prey and predatorboth present stably oscillatory motions with respect totime, which is equivalent to closed orbits in the phasespace. Using this deterministic approach, there is a sta-bility indicator χ corresponding to each orbit(Fig.(3)).

χ = −by − dx+ a log y + c log x (4)

However, the replicator equations couldn’t capture whattruly happens for a finite population system. There arealways intrinsic fluctuations, which can not be neglected,due to finite population. What is the outcome of thisstochastic behavior for Lotka-Volterra model? In Fig.(3),we observe closed circles with positive curvature in phasespace. Therefore, from our method, when noise is pre-sented in this situation, dissipative motion would appear.The system will deviate away from its original contourand finally go to extinction because of absorbing bound-aries in phase space. In order to see how the methodagrees with realistic situation, we calculate the first mo-ment of χ, which indicates dissipation, in phase space.The algorism for the stochastic motion is to use proba-bility to describe the mechanism. Not surprisingly, χ isnegative over the whole phase space. Fig.(4)shows greatagreement between the algorism and the theoretical cal-culation in region A of Fig.(3). In region A, it’s nat-ural to have successful agreement because that all theapproximations we have done during calculation are sat-isfied. During the derivation of the method, there are two

assumptions we have used. One is that the standard de-viation of Gaussian distribution of noise is much smallerthan the radius of curvature of contour. The other is thatduring the hopping of species from one point to anotherin the phase space, gradient of stability indicator and cur-vature of contour remain the same within the range of thestandard deviation. In region A, the curvature of contouris small compared to other areas and the change of cur-vature and of gradient is also small. Therefor, a gooddemonstration between algorism and theoretical calcula-tion is presented in this area. However, the conditions inregion B in Fig.(3) are far away from the approximationsand hence conflictions appear.

The neutral coexistence predicted by deterministic ap-proach is broken due to the intrinsic fluctuation of a bi-ological system. The fluctuation-induced dissipation isgenerated and we use the first moment of χ to capturethe dissipative behavior. Successfully, the great agree-ment between the algorism and the theoretical calcula-tion in the area we focus on is found in Lotka-Volterramodel and this presents the validity of the method.summary. We present a method that can concretelypredict the fluctuation-induced dissipation in a biologi-cal system. Although the deterministic approach illus-trates the beautiful coexistence, this method ignores theunavoidable noise in the evolution due to the externalinfluence of environment or the intrinsic fluctuation ofa finite population system. As a result, it is better touse stochastic method to describe evolution; however,the extinction of the species always occur in this kindof simulation. To have a more profound understandingof the extinction process, we take an appropriate coarse-graining of phase space structure and define a stabilityindicator χ to describe the extinction process of system,finding that the fluctuations always induce the dissipa-tion by the positive curvature of the equi-χ contour ei-ther in neutral stable or stable cases. As a demonstratingexample, we compare theoretical calculation and simula-tion of stochastic Lotka-Volterra model and find fantasticagreement. Furthermore, our finding shares great simi-larity with fluctuation-dissipation theorem in statisticalphysics while the major difference is that our method isvalid even in the system which is far from equilibrium.The breakdown of biological diversity in ecosystem is notso mysterious. From the geometric structure of the sta-bility indicator χ in phase space, we can see that thefluctuation-induced dissipation is predicted and plays animportant role in the evolutionary dynamics for biologi-cal systems.

DISCUSSION

The spirit of the our idea is that fluctuation can in-duce a non-zero dissipation through the geometric struc-ture of the phase space determined by the determinis-

4

tic dynamics. Therefore, the idea is close related to thefluctuation-dissipation theorem in statistical physics. Tofurther illustrate the similarity, for two-species system,under the the same approximation as in the calculationof the first moment, we can derive the formula for thesecond moment of stochastic variable χ:⟨

∆χ2⟩

= |∇χ|2σ2 (5)

which is a reasonable result because second moment re-flects magnitude of fluctuation of the stochastic variableχ which depends on the density of contours and inten-sity of fluctuation only, instead of curvature of phase con-tours.

Take the continuous limit of time interval, the correla-tion of noise on x1 and x2 is 〈ξi(t1)ξi(t2)〉 = Dδ(t1 − t2),i = 1, 2, then we can derive the generalized Langevinequation for the stochastic variable χ:

dχ

dt= −1

2κ|∇χ|D + ξχ(t) (6)

where 〈ξχ(t1)ξχ(t2)〉 = |∇χ|2Dδ(t1 − t2). From this ex-pression, it is obvious to see that the dissipative dynamicsemerges from the seemingly harmless intrinsic fluctua-tions, and the connection of fluctuation and dissipationis in the same way as Fluctuation-dissipation theorem.However, the marked difference is that the relation wederive can apply to the system which is far from equilib-rium while Fluctuation-dissipation theorem only worksin the linear response regime of a system. Since theidea is valid for the system which is far from equilib-rium, we may interpret our result as a generalization ofFluctuation-dissipation theorem for the ecosystem.

In the above discussion, we focus on a system where allspecies display dynamical balance with a neutral stabil-ity, and discover that the fluctuations induce a dissipa-tion which drives species to extinction. However, most ofdynamical systems show either stable or unstable stabil-ity instead of a neutral dynamical balance. While extinc-tion in the unstable system is expected, the extinctionprocess of a stable system is more complicated becauseit is hard to capture the effect caused by fluctuations. Inorder to investigate the effect of fluctuation on the sta-ble system, consider a N-species system where there isonly one stable fixed point. Although the deterministicflow in the phase space spanned by population numbersof species xi(i = 1, 2...N) drives the system toward thefixed point, we can naively think that extinction alwaysoccurs due to a series of unfortunate events as long asthe evolution time is long enough. To measure the ef-fect caused by fluctuations quantitatively, we can definestability indicator χ as the following:

χ(x1, x2, ...xN ) = −αN∑i=1

(xi − x∗i )2 (7)

where (x∗1, x∗2, ...x

∗N ) denotes the coordinate of fixed

point, and α is a positive constant. By this defini-tion, the set of equi-χ points forms a closed hypersurfacewhich centers at the fixed point, where χ achieves maxi-mum. Therefore, we can measure the extinction processby studying the dissipation of χ, which indicates the evo-lution of the system toward the boundaries. Ignore thenoise first, different from the case of neutral dynamicalsystem, the deterministic flows contribute to the timeevolution of χ as well. For the effect of fluctuations, be-cause the equi-χ hypersurface is closed, we can still applyour previous method to determine the dissipation causedby fluctuations. Again, during a small time interval, theset of points with the same probabilities that system canhop on forms a closed hypersurface. With the positivecurvature of equi-χ hypersurface, there is an asymmetryof phase space volume between inner and outer region,and thus the fluctuation-induced dissipation is expected.Furthermore, we can derive the Langevin equation of χ:

dχ

dt= ∇χ · (d~x

dt)flow −

1

2κ|∇χ|D + ξχ(t) (8)

where 〈ξχ(t1)ξχ(t2)〉 = |∇χ|2Dδ(t1 − t2), and (d~xdt )flowdenotes the deterministic flow in the phase space whichalways increases the value of χ. Therefore, the extinctionprocess of species is complicated due to the competitionbetween deterministic flow and fluctuation-induced dis-sipation. However, due to the existence of white noise,the probability distribution in the χ will spread through-out the χ space. In consequence, the system will still goto extinction due to the absorbing boundaries of phasespace spanned by N-species for a large enough time.

METHODS

Differential equations in continuous limit. Themost fundamental description for biological evolution isto depict the discrete change of species population withprobabilities. From this microscopic view of how thespecies evolves, Eq.(1) can be derived. To start with,we take the Lotka-Volterra model, with two species X1

and X2, as an example. For species X1, during the sim-ulation time interval ∆τ , where τ = 1, 2, 3, ..., the popu-lation change ∆x1 may be 1, 0 or −1, which correspondsto the probability P+1, P0 and P−1 respectively. We canwrite down the average change of x1 during a simulationtime interval.

〈∆x1〉 = 1 · P+1 + 0 · P0 + (−1) · P−1= ax1∆t− bx2x1∆t (9)

Probability can be naturally defined as Pα = Aαx1∆t,where Aα depends on evolution condition. For example,P+1 = a and P−1 = bx2. Besides, to make the valueof probability smaller than one, we set the natural time

5

FIG. 3: Geometric definitions for the angular variables usedin the derivations.

interval ∆t = ∆τ/N2 = 1/N2 in this case. Since Eq.(18)is the average value, to have the explicit expression of∆x1 during ∆τ , the noise must be added.

∆x1 = 〈∆x1〉+ η1 (10)

Next, we devide by ∆τ and use the property ∆τ = 1:

∆x1∆τ

=〈∆x1〉

∆τ+

η1∆τ

= 〈∆x1〉+ η1 (11)

Using the relation ∆τ = N2∆t, we replace ∆τ in theleft hand side with N2∆t and move N2 to the right.Therefore, we have the following equation:

dx1dt' ∆x1

∆t= 〈∆x1〉N2 + η1N

2

= ax1 − bx2x1 + ξ1 (12)

When populatin number N is large enough, ∆t would belimited to zero and hence we can approximate the discretepopulation change to continuous differential equation. IfEq.(21)is taken the average, the noise term would disap-pear and then turn back to the replicator equation. Thesimilar derivation could be applied to multiple species.As a result, we have the generalized version Eq.(1)in thefirst place. This part of demonstration shows that thestochastic differential Eq.(1) is derived from the micro-scopic view of how species interacts with others. Withlarge population N approximation, we have the differ-ential equation which comes from the discrete stochasticequation.Dissipation in biodiversity indicator. Here wepresent the details for the calculation of 〈∆χ〉 and

⟨∆χ2

⟩:

For system located at point A as shown in (Fig.3), if wemark all the probable points that system can hop withinsmall time interval ∆t, the set of points will form a cir-cle centered at point A under the approximation thatrandom white noises obey Gaussian distribution with

equal standard deviation. Inside the circle with radiusr, assuming |∇χ| and κ are constant, we can calcu-late the change of χ corresponding to different hoppingpoints∆χ(r, θ):

∆χ(r, θ) = −|∇χ|AC = −|∇χ|(R−R cosφ+ r cos θ)(13)

Use the distance between point B and AC to find therelation between φ and θ

R sinφ = r sin θ

cosφ =

√R2 − r2 sin2 θ

R

(14)

Substituting Eq.(14)to Eq.(13):

∆χ = −|∇χ|(R−√R2 − r2 sin2 θ + r cos θ) (15)

Taking the angular average of ∆χ:

〈∆χ〉θ =1

2π

∫ 2π

0

∆χdθ

=1

2π

∫ 2π

0

−|∇χ|(R−√R2 − r2 sin2 θ + r cos θ)dθ

= −|∇χ|2π

[2πR−R

∫ 2π

0

√1− (

r sin θ

R)2dθ

](16)

under the approximation that r � R, we can evaluatethe integral by binomial expansion to first order:∫ 2π

0

√1− (

r sin θ

R)2dθ =

∫ 2π

0

(1− r2

2R2sin2 θ

)dθ

= 2π − |∇χ|r2

4R(17)

Substituting Eq.(17) to Eq.(16):

〈∆χ〉θ = −|∇χ|r2

4R= −1

4κ|∇χ|r2 (18)

Following similar steps, procedure, we can calculate

⟨∆χ2

⟩θ

= 〈∆χ〉θ =1

2π

∫ 2π

0

dθ∆χ2 =1

2|∇χ|2r2 (19)

From Eq.(18) and Eq.(19), we can already observe thatthe first moment has linear dependence on curvature andgradient χ while the second moment depends on gradientonly. To obtain explicit result, assuming the randomdisplacement(x, y) caused by random white noises obeyGaussian distribution,

Px(x) =1√

2πσ2e−

x2

2σ2

Py(y) =1√

2πσ2e−

y2

2σ2 .

(20)

6

It is straightforward to compute the radial distributionof noises,

Pr(r) =r

σ2e−

r2

2σ2 (21)

where S refers to the area in the circle with radius r.Assuming σ � R, due to extreme narrow peak of prob-ability distribution, we can safely use Eq.(18) in the cal-culation of radial average of 〈∆χ〉θ even for r → ∞ inthe integral . Thus we can evaluate ∆χ by taking ra-dial average of 〈∆χ〉θ from the contributions of differenthopping radius:

〈∆χ〉 =

∫ ∞0

drPr(r) 〈∆χ〉θ

= − κ

4σ2|∇χ|

∫ ∞0

r3e−r2

2σ2 dr

= −1

2κσ2|∇χ|

(22)

Similarly, the average of the second moment can also beobtained, ⟨

∆χ2⟩

= |∇χ|2σ2. (23)

[1] R. M. May, Stability and Complexity in Model Ecosys-tems (Princeton University Press, Princeton, New Jersey,1974), second edition.

[2] S. L. Pimm, Nature 307, 321 (1984).[3] D. Jablonski Proc. Natl. Acad. Sci. USA 105, 11528

(2008).[4] J. F. McLaughlin, J. J. Hellmann, C. L. Boggs, and P.

R. Ehrlich, Proc. Natl Acad. Sci. USA 99, 6070 (2002).[5] C. Both, S. Bouwhuis, C. M. Lessells, and M. E. Visser,

Nature 441, 81 (2006).[6] O. E. Sala et al., Science 287, 1770 (2000)[7] T. Reichenbach, M. Mobilia and E. Frey, Nature 448,

1046 (2007).[8] M. Loreau et al., Science 294, 804 (2001).[9] J. M. Smith, Evolution and the Theory of Games (Cam-

bridge Univ. Press, Cambridge, 1982).[10] J. Hofbauer and K. Sigmund, Evolutionary Games and

Population Dynamics (Cambridge Univ. Press, Cam-bridge, 1998).

[11] M. A. Nowak, Evolutionary Dynamics (Belknap Press,Cambridge, Massachusetts, 2006).

[12] M. A. Nowak, Science 314, 1560 (2006).[13] B. Drossel, Adv. Phys. 50, 209 (2001).[14] B. Kerr, M. A. Riley, M. W. Feldma and B. J. M. Bo-

hannan, Nature 418, 171 (2002).[15] T. L. Czaran, R. F. Hoekstra and L. Pagie, Proc. Natl

Acad. Sci. USA 99, 786 (2002).[16] M. A. Nowak and K. Sigmund, Science 303, 793 (2004).[17] S. A. West, A. S. Griffin, A. Gardner and S. P. Diggle,

Nature Rev. Micro. 4, 597 (2006).[18] G. Szabo, and G. Fath, Phys. Rep. 446, 97 (2007).[19] A. Traulsen, J. C. Claussen and C. Hauert, Phys. Rev.

Lett. 95, 238701 (2005).

[20] T. Reichenbach and E. Frey, Phys. Rev. Lett. 101,058102 (2008).

[21] J. C. Claussen and A. Traulsen, Phys. Rev. Lett. 100,058104 (2008).

[22] M. Berr, T. Reichenbach, M. Schottenloher and E. Frey,Phys. Rev. Lett. 102, 048102 (2009).

[23] J. B. C. Jackson and L. Buss, Proc. Natl Acad. Sci. USA72, 5160 (1975).

[24] B. Sinervo and C. M. Lively, Nature 380, 240 (1996).[25] B. C. Kirkup and M. A. Riley, Nature 428, 412 (2004).[26] M. Ifti and B. Bergersen, Eur. Phys. J. E 10, 241 (2003);

Eur. Phys. J. B 37, 101 (2004).[27] T. Reichenbach, M. Mobilia, and E. Frey, Phys. Rev. E

74, 051907 (2006).

ACKNOWLEDGEMENTS

We acknowledge supports from the Ministry of Scienceand Technology in Taiwan through grant MOST 103-2112-M-007-011-MY3. Financial supports and friendlyenvironment provided by the National Center for Theo-retical Sciences in Taiwan are also greatly appreciated.

AUTHOR CONTRIBUTIONS

T.C.L. and Y.K.C. performed the analytical and nu-merical calculations. H.H.L. and C.C.C. supervise thewhole work. All authors contributed to the preparationof this manuscript.