Simulation and Calibration Between Parameters of Continuous Time Random Walks and Subdifusive Model ^†

ABSTRACT

We address the problem of subdiffusion or normal diffusion to perform a calibration between simulations parameters and those from a subdiffusive model. The theoretical model consists to a generalized diffusion equation with fractional derivatives in time. The data generated by simulations represents continuous-time random walks with controlled mean waiting time and jump length variance to provide a full range of cases between subdiffusion and normal diffusion. From simulations, we compare the accuracy of two methods to obtain the diffusion constant and the order of fractional derivatives: the analysis of the dispersion of the variance in time and an optimized fitting of the histograms of positions with theoretical model solutions. We highlight the connection between the parameters of the simulations those of theoretical models.

Keywords:
anomalous diffusion; fractional diffusion equation; calibration

RESUMO

Abordamos o problema da subdifusão e da difusão normal a fim de realizar uma calibração entre parâmetros simulados e parâmetros de um modelo subdifusivo. O modelo teórico consiste em uma equação de difusão generalizada por derivadas fracionárias temporais. Os dados são gerados por simulações baseadas em caminhadas aleatórias de tempo contínuo com média do tempo de espera e variância do comprimento de saltos controlados para fornecer uma gama completa de casos entre a subdifusão e a difusão normal. A partir das simulações comparamos a precisão entre dois métodos para obter os parâmetros coeficiente de difusão e ordem da derivada fracionária: análise da dispersão da variância no tempo e ajuste por otimização dos histogramas de posições com as soluções do modelo teórico. Destacamos a conexão entre os parâmetros das simulações e dos parâmetros do modelo teórico.

Palavras-chave:
difusão anômala; equação de difusão fracionária; calibração

1 INTRODUCTION

Anomalous diffusion has been attracting the attention of the scientific community due the large plethora of natural systems which display large deviations from normal diffusion, as plasma diffusion ³, fluid flow in porous media ^{3), (22}, diffusion in fractal structures ²³, turbulence ⁸ etc. Differently to the normal diffusion, associated to local, short range correlations, the presence of long range correlations lead to non-Gaussian probability density functions (PDF).

In the class of problems generically called anomalous diffusionsome or anomalous Brownian motions, the microscopic kinetics is characterized by increasing probabilities to observe large steps, causing non-Gaussian PDF with heavy tales, or conversely, to observe large waiting times elapsed between steps causing non-Gaussian PDF, causing non-Gaussian PDF with light tales. From a macroscopic point-of-view, those anomalous diffusion processes do not fit in the classical diffusion description by a Fokker-Planck equation for the evolution of the PDF u(x, t) in space and time ^{11 24}.

As a direct consequence of the central limit theorem, Markov processes displays a linear evolution of the mean square displacement with time, 〈x ²(t)〉 ∝ Dt, featuring the normal, or Brownian, diffusion ²⁵. According to that theorem, if the single variable PDF has finite second moment, the cumulative density functions (CDF) is also gaussian.

On the other hand, if the PDF has long-range tails, determining the divergence of the second moment, the Lévy-Gnedenko limit theorem states that the CDF is Lévy distributed, meaning that it has long-range tails too ⁹. Alternatively, the q-central limit theorem states that, in the same conditions for the PDF, a q-Gaussian CDF with 1 < q < 3 is obtained ²⁷. Those process are generally characterized by a non-linear dependence of the mean square displacement in time, with 〈x ²(t)〉 ∝ Dt ^α17. According to the value of the exponent α, the anomalous processes can be classified as sub-diffusive (0 < α < 1) or super-diffusive (1 < α < 2) ⁽¹⁷. A relation between the diffusion exponent α and the q parameter was proposed by Bukman and Tsallis ²⁶, and recently verified experimentally ⁶.

Random walks flight simulations can be generalized as continuous time random walk (CTRW) which are able to simulate processes displaying anomalous diffusion. In particular, considering a random distribution of step lengths as well a random distribution on waiting times, it is possible to control the process in order to simulate both sub-diffusion or super-diffusion regimes ^17),18),19.

The description of the macroscopic behavior of u(x, t) for anomalous diffusion process can be modeled by means of generalized diffusion equations with fractional derivatives in time as well as in space ^{1),2),13),14),15}. Another class of generalizations of diffusion equations is obtained by the inclusion of nonlinearity in the partial differential equation ^4),13),25. Those complementary approaches have motivated this study, which aims understanding the phenomena of anomalous diffusion in complex dynamic systems ^{14),16),18),20}.

From one side, there is the general interest in simulating microscopic anomalous diffusive systems in order to reproduce some observable macroscopic behavior with its associated variables. By another side, there lies the interest to constitute macroscopic anomalous diffusive models to describe some microscopic systems with detailed dynamics not easily observable. Thus, a question arises naturally: is there a systematic relation between microscopic simulation parameters and the macroscopic anomalous diffusion model? In the case of a positive answer, which is that relation? Moreover, one can wonder which is the best method to compare models to simulations as well as which models are more appropriate to fit simulations.

In the present work, we address the problem of subdiffusion and normal diffusion to perform a calibration between simulation microscopic parameters and macroscopic theoretical model parameters of a generalized anomalous diffusion (subdiffusion) model. The data is generated by simulations of CTRW with controlled mean waiting time and jump length variance to provide a full range of cases between subdiffusion and normal diffusion. The theoretical model consists in a contour value problem with asymptomatically vanishing boundary conditions having a fractional partial differential equation (FPDE) as an subdiffusion equation, with fractional derivatives in time. From the simulations, we obtain the diffusion constant and the order of fractional derivatives by two different approaches: (1) optimization fitting of theoretical model solutions to the histogram of positions of the particles and (2) by analysis of the dispersion in time of the variance of the particles position. The microscopic parameters of the simulations are then associated to the theoretical macroscopic parameters. We then highlight the connection between the parameters of the simulations the parameters of the theoretical models for the cases of fractional derivatives in time. In Section (2), we present the generalized model and solutions. The long range behavior of the solutions is discussed in Section (3). Section (4) comprises the aspects of the simulations which is followed by the Section (5) of results and Section (6) for the conclusions.

2 SUBDIFFUSION EQUATION AND INITIAL BOUNDARY VALUE PROBLEM

To theoretically describe the subdiffusion, the differential partial equation for the probability density function u(x, t) of finding a particle between x and x + δx at the time t is written by means of a generalized diffusion equation with fractional derivatives in time. Hence, a uni-dimensional open contour value problem subject to asymptotic contour conditions and appropriated initial conditions can be written as

$\begin{array}{l}{}^C\mathcal{D}_t^{\alpha }u(x,t)=D_{\alpha }\frac{{\partial }^{2}u(x,t)}{\partial x^2},\\ u(\pm \infty ,t)=0,u(x,0)=N\delta (x-x_0),u_t(x,0)=0,\\ -\infty <x<\infty t>0\end{array}$ (2.1)

Where _^Dα is the subdiffusion coefficient corresponding to constant of proportionality between the concavity of the solution and the generalized velocity of the process represented by the fractional derivative with respect to time. In turn α represents the order of this derivative. The fractional derivative operator in time ${}^C\mathcal{D}_t^{\alpha }$ is defined in the Caputo sense in order to allow the statement of initial conditions with physically meaningful derivatives of integer order in time ²¹,

${}^C\mathcal{D}_t^{\alpha }f=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t{f}^{\text{'}}\left(\tau \right)\frac{1}{(t-\tau {)}^{\alpha }}d\tau ,0<\alpha <1.$ (2.2)

For α = 1 in (2.1), the classical diffusion equation is recovered. With varying the fractional order of the derivatives, 0 < α < 1, one can describe subdiffusive processes ¹⁷.

The solution of the problem (2.1) reads

$u(x,t)=\frac{N}{2\sqrt{D_{\alpha }}{t}^{\frac{\alpha }{2}}}\sum _{k=0}^{\infty }\frac{(-z{)}^k}{k!\Gamma (-\frac{\alpha }{2}k+(1-\frac{\alpha }{2}\left)\right)},z=\frac{|x-x_0|}{\sqrt{D_{\alpha }}{t}^{\frac{\alpha }{2}}}.$ (2.3)

Such solution describe subdiffusive process. After normalization of the solution (2.3) one finds the Mittag-Leffler PDF _^Gα (x, t) ¹⁰. By using (2.1) and (2.2), one can obtain the following expression for the mean square displacement:

$⟨x^2\left(t\right){⟩}_G=\frac{2D_{\alpha }}{\Gamma (\alpha +1)}{t}^{\alpha },0<\alpha <1.$ (2.4)

We use these case to fit the histogram of displacements of CTRW simulations to produce subdiffusive processes.

3 LONG RANGE BEHAVIOR OF THE MODEL

By means of the the Laplace transform in the Caputo sense on the time variable ℒ{u(x, t); s} = ũ (x, s) and the Fourier transform on space variable, ℱ{u(x, t); k}=û (k, t), one obtains from (2.1) its Laplace-Fourier transform ^{18), (21},

$\hat{\tilde{u}}(k,s)=\frac{s^{\alpha }\hat{u}(k,0)}{s(s^{\alpha }+D_{\alpha }{k}^2)}.$ (3.1)

For a CTRW, the probability distribution function of a walk to the a given position x in a time t is given by the length of jump as well as the waiting time between successive jumps are well described by a function Ψ(x, t), whose the marginal probabilities in the space and time are respectively given by λ(x) and ψ(t) ¹⁷. It is possible to show that the Fourier-Laplace transform of the PDF u(x, t) of finding the particle in the time t , between x and x + dx is given by

$\hat{\tilde{u}}(k,s)=\frac{1-\tilde{\psi }\left(s\right)}{s}\frac{\tilde{\hat{u}}(k,0)}{(1-\tilde{\psi }(s\left)\hat{\lambda }\right(k\left)\right)},$ (3.2)

which is known as the Montroll-Weiss formula ^12),19. The long range behaviour of above equation is accessed in the infrared limit of the Fourier-Laplace transform, s → 0, k → 0. By representing the caracteristic functions, τ and σ ², in moments expansion by using, respectively, Fourier and Laplace transform of λ(x) and ψ(t) one finds

$\begin{array}{l}\tilde{\psi }\left(s\right)\approx 1-{\tau }^{\alpha }{s}^{\alpha },0<\alpha \le 1,\\ \hat{\lambda }\left(k\right)\approx 1-{\sigma }^2{k}^2\end{array}$ (3.3)

where, τ and σ are scale constants such that τ ^α and σ ² are identified to the first and second moments when α = 1 ¹⁷. The equations (3.3) describes the Fourier-Laplace transform of the PDF in which the first moment, τ, diverges while the second moment, σ ², converges. To the first order in (3.3), equation (3.2) reduces to

$\hat{\tilde{u}}(k,s)=\frac{s^{\alpha }\tilde{\hat{u}}(k,0)}{s(s^{\alpha }+\frac{{\sigma }^2}{2{\tau }^{\alpha }}{k}^2+\cdots )}.$ (3.4)

The comparison between _^Dα of equation (3.1) and the scales σ and τ above allows us to obtain the equation

$\log \left(2D_{\alpha }\right)=\log \left({\sigma }^2\right)-\alpha \log \left(\tau \right).$ (3.5)

These equation exhibits a scale relation between the moments arising from the simulation and macroscopic parameters of the theoretical model.

Let the distributions of jump length δx and waiting times δt be, respectively, λ'(δx) and ψ'(δt). Then, the type of anomalous diffusion simulated by CTRW is also determined also by the long range behavior of first moment 〈δt〉_ψ' and the second moment 〈(δx)²〉_λ' by the following rule ¹⁹:

• Case 1. If 〈δt〉_ψ' →∞ and 〈(δx)²〉_λ' is finite, the process is subdifusive.

• Case 2. If 〈δt〉_ψ' is finite and 〈(δx)²〉_λ' →∞, the process is superdifusive.

• Case 3. If 〈δt〉_ψ' and 〈(δx)²〉_λ' are both finite, the process is a normal diffusion.

Such cases can be read as prescriptions to produce CTRWs with the desired types of anomalous diffusions. The case 2, superdiffusive, presents several subtleties and characteristics that would lengthen the discussion and escape the scope of the present manuscript. In this way we decided to leave this analysis for future work and restrict ourselves to the cases 1 and 3, subdiffusive and normal diffusion.

4 SIMULATION OF ANOMALOUS RANDOM WALKS

A convenient distribution to control divergence of the moments is the power law PDF, K (ξ) = $\frac{a}{(\xi +\epsilon {)}^p}$, having a as a normalization constant, ξ ≥ 0, p > 1, Λ stands for the upper limit cutoff for the average integrals, and ε, the lower limit one. In considering jump length, δx, and waiting time, δt, PDFs separated, we have

${\lambda }^{\text{'}}\left(\delta x\right)=\frac{A}{\left(\right|\delta x|+\epsilon {)}^r},and{\psi }^{\text{'}}\left(\delta t\right)=\frac{B}{(\delta t+\epsilon {)}^s}.$ (4.1)

After determining a, the m-th moment reads

$⟨{\xi }^m⟩=\frac{p-1}{p-(m+1)}{\epsilon }^m\left[\frac{1-{\left(\frac{\epsilon }{\Lambda +\epsilon }\right)}^{p-(m+1)}}{1-{\left(\frac{\epsilon }{\Lambda +\epsilon }\right)}^{p-1}}\right].$ (4.2)

Formally, if p > m + 1, 〈ξ^m 〉 is convergent in the limit Λ → ∞. By setting the power dependencies r and s with the criterion of the divergence p > m + 1 from equation (4.2), and consequently, determining the convergence or divergence of the first moment 〈δt〉_ψ' or of the second moment 〈(δx)²〉_λ', it is possible to fulfill the cases 1 to 3 in Section (3), controlling the anomalous diffusion processes ranging from subdiffusion to superdiffusion.

Simulated processes consisted of unidimensional CTRW composed by a population of N = 10.000 particles with step lengths and waiting times randomly determined by a Monte Carlo algorithm from power law PDFs (4.1). The time stop criterion was chosen to be so that Σ δt ≥ 2¹⁰, since the number of simulations showed that the parameters remained almost constant. Figure (1) shows the time evolution of the population variance, log₂(t)×log₂(〈x ²(t)〉) for cases of subdiffusion and normal diffusion. The histogram of positions of the particles is depicted in Figure (2) for the subdiffusion case.

Figure 1:
Time evolution of log₂(〈x ²(t)〉) generated by CTRW of a population of N = 10.000 particles over a time t ≥ 2¹⁰ with random jump lengths and random waiting times, both obtained by power law PDF with r = 4 and different values powers s. According to the parameters, the process is set to be normal or subdiffusive. The lines are fitted by linear regression in order to determine the model parameters. Case s = 3 and r = 4 correspond to a normal diffusion process.

Figure 2:
Histogram of positions of a population of of N=10.000 particles generated by simulation over a time t ≥ 2¹⁰ of a CTRW with random jump lengths and random waiting times, both obtained by power law PDF. According to the simulation parameters r and s, the process is set to be normal or subdiffusive.

5 RESULTS

Simulations were conducted by defining the parameters of the power laws PDF for jump length and waiting times. According to the definition of the parameters and, as the aforementioned cases 1 and 3 (Section 3), the processes were set to be subdiffusive or normal. The assessment of the macroscopic parameters from the simulated microscopic data which was acquired by two different approaches. Firstly, we analyzed the dispersion of the particles population by applying linear regression over the data in the space log₂(t)×log₂(〈x ²(t)〉) according to the equation (2.4). The parameters α and ^_Dα are achieved from subdiffusive processes simulations and are given, respectively, by the slope and ordinate intercept of the linear regression plots between time, t, and variance, 〈x ²〉, (see Fig. 1). Since those parameters are obtained by the simulated data without any reference to the theoretical model, they receive the subscript S assigning for ‘simulation’. The results of analysis of studied cases can be seen in Table (5).

Table 1:
Simulation data generated by power law PDFs for random jump sizes with r = 4 and for waiting time varying s in order to provide subdiffusion processes and the respective values of the parameters α and D _α . The parameters assigned by ’S’ are obtained by linear regression over the time evolution of the variance and those assigned by ’T’ are obtained by fitting the model solution G _α (x, t) (2.3) by BFGS optimization technique. The mean squared deviations MSDs are also calculated in both cases.

s	αS	$D_{{\alpha }_S}$	αT	$D_{{\alpha }_T}$	MDSS	MDST
3.0	0,99212726	0,42963683	0,99999965	0,41181830	10-9	10-9
1.9	0,84902374	0,63356256	0,84358289	0,63223199	10-8	10-8
1.8	0,78971432	0,72844909	0,78471521	0,72715277	10-8	10-8
1.7	0,70974939	0,94427836	0,69789268	0,94184530	10-8	10-8
1.6	0,60320980	1,34875574	0,58446281	1,34542329	10-7	10-7
1.5	0,48760105	1,91351808	0,47772984	1,91212494	10-6	10-6
1.4	0,38735396	2,52506216	0,38435808	2,52369939	10-6	10-6
1.3	0,31066169	3,22746669	0,31661584	3,22746669	10-7	10-7
1.2	0,25982419	3,76553682	0,27233279	3,76625088	10-7	10-7
1.1	0,23832064	4,28744069	0,25122937	4,27971155	10-6	10-6

The second method to obtain the macroscopic parameters comprises the optimization fitting of the model solution _^Gα (x, t) (2.3) to the histogram of position of the processes according to the subdiffusion. Hence, _^Gα (x, t) is fitted to subdiffusion processes to asses the parameters α and _^Dα . The adopted optimization method was the Broyden-Fletcher-Goldfarb-Shanno (BFGS) ⁵ using a minimization function that considers the mean square error of the deviations. The parameters obtained by that approach receive the subscript T in mention to the theoretical solutions. Results of such analysis are depicted in Figure (3) and the parameters and mean squared deviations (MSD) are listed in Table (5).

Figure 3:
Model solutions $G_{{\alpha }_T}$(x, t) (2.3) with parameters assessed by optimization fitting to the histogram of positions generated by simulation over a time t ≥ 2¹⁰ of a CTRW for a population of N = 10.000 particles.

From the MSD in the Table (5), one can notice that the methods of optimization fit of the histogram and the method of linear regression from the time evolution of the variance were similar in accuracy. The simulations have produced processes with the expected behavior as subdiffusion or normal diffusion in accord to the prescriptions given in cases 1 and 3 (Section 3) and to the criterion p < m + 1(r < 3, s < 2).

A calibration on the parameter space between simulation and theoretical model is depicted in (4) for varying population N. The line slope is compatible with the effective values of τ* = -4.4935 and σ ² = 4.1064 in equation (2.4).

Figure 4:
Calibration between microscopic or simulation parameters and macroscopic or theoretical model parameters.

6 CONCLUSION

Both methods, in order to determine the parameters of the model, presented the same mean squared deviations. The method of the fitting dispersion by the time evolution of the variance, however, produces results with much lower computational cost and presents the same dispersion of the optimized fitting histograms of positions.

The definition of the power laws is not uniquely dependent of the power parameter p, (r, s), but it also demands the choose of two additional cutoffs; the first ε to avoid the singularity at zero and the upper limit cutoff Λ to provide integrability when a given moment of interest (σ ², τ) is divergent. In order to avoid excessive arbitrariness, we make a prescription which is to set a constraint between both cutoffs. Under such condition, the divergent factor of the moment of interest becomes near the unit.

It was possible to determine a set of equations and prescriptions to estimate which are the powers of the power law distributions for the jumps and waiting time would produce the macroscopic parameters of the linear model written in terms of partial differential equations with fractional derivatives. However, not all cases of fractional derivatives in time were explored. We restricted ourselves to the case which have one of the derivatives being of integer order: the second order derivative in space was fixed and the time derivative remained fractional. Although we could generate a table with those values. We found a strong correlation between the values of the simulation parameters and those of the theoretical parameters.

REFERENCES

• ¹
  A.A.M. Arafa & S.Z. Rida. Exact solutions of fractional-order biological. Communications in Theoretical Physics, 6 (2009), 992-996.
• ²
  O.P. Argrawal. Fractional variational calculus in terms of riesz fractional derivatives. J. Phys. A: Math. Theor., 40 (2007), 6287-6303.
• ³
  J.G. Berryman. Evolution of a stable profile for a class of nonlinear diffusion equations with fixed boundaries. Journal of mathematical physics, 18 (1977), 2108-2115.
• ⁴
  M. Bologna, C. Tsallis & P. Grigolini. Anomalous diffusion associated with nonlinear fractional derivative fokker-planck-like equation: exact time-dependent solutions. Physical Review E, 62 (2000), 2213.
• ⁵
  R.H. Byrd, P. Lu & J. Nocedal. A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing, 16 (1995), 1190-1208.
• ⁶
  G. Combe, V. Richefeu, M. Stasiak & A.P.F. Atman. Experimental validation of a nonextensive scaling law in confined granular media. Phys. Rev. Lett., 115 (2015), 238301.
• ⁷
  G.L. Eyink & H. Spohn. Negative-temperature states and large-scale, long-lived vortices in two-dimensional turbulence. Journal of statistical physics, 70 (1993), 833-886.
• ⁸
  G.A. Garosi, G. Bekefi & M. Schulz. Anomalous diffusion and resistivity of a turbulence, weakly ionized plasma. Appl. Phys. Lett., 15 (1969), 334-337.
• ⁹
  B.V. Gnedenko & A.N. Kolmogorov. Limit distributions for sums of independent random variables. Bull. Amer. Math. Soc., 62 (1956), 50-52.
• ¹⁰
  H.J. Haubold, A.M. Mathai & R.K. Saxena. Mittag-leffler functions and their applications. Journal of Applied Mathematics, 2011 (2011), 51 pp.
• ¹¹
  R. Hilfer, editor. Applications of fractional calculus in physics. World Scientific, (2000).
• ¹²
  V.M. Kenkre, E.W. Montroll & M.F. Shlesinger. Generalized master equations for continuous-times random walks. Journal of Statistical Physics, 9 (1973), 45-50.
• ¹³
  E.K. Lenzi, R.S. Mendes & C. Tsallis. Crossover in diffusion equation: Anomalous and normal behaviors. Physical Review E, 67(3) (2003), p. 031104.
• ¹⁴
  F. Mainardi. Fractional relaxation-oscillation and fractional diffusion-wave phenomena. Chaos Solitons & Fractals, 7 (1996), 1461-1477.
• ¹⁵
  F. Mainardi. The fundamental solutions for the fractional diffusion-wave equations. Applied Mathematics Letters, 9(6) (1996), 23-28.
• ¹⁶
  F. Mainardi. Fractional Calculus and Waves in Linear Viscoelasticity: An Introduction to Mathematical Models. Imperial College Press, 2nd edition, (2010).
• ¹⁷
  R. Metzler & J. Klafter. The random walk’s guide to anomalous diffusion: A fractional dynamics approach. Physics Reports, 339 (2000), 1-77.
• ¹⁸
  R. Metzler & T.F. Nonnenmacher. Space- and time-fracitonal diffusion and wave equations, fractional fokker-planck equations, and physical motivation. Chemical Physics, 284 (2002), 67-90.
• ¹⁹
  E.W. Montroll & G.H. Weiss. Random walks on lattices II. J. Math. Phys., 6 (1965), 167.
• ²⁰
  Paolo Paradisi, Rita Cesari, F. Mainardi & Francesco Tampieri. The fractional fick’s law for non-local transport processes. Physica A: Statistical Mechanics and its Applications, 293 (2001), 130-142.
• ²¹
  I. Podlubny. Fractional Differential Equations. Academic Press, (1999).
• ²²
  M.F. Shlesinger & J. Klafter. On the relationship among three theories of relaxation in disordered systems. Proceedings of the National Academy of Sciences, 83 (1986), 848-851.
• ²³
  J. Stephenson. Some non-linear diffusion equations and fractal diffusion. Physica A, 222 (1995), 234-247.
• ²⁴
  L.T. Takahashi et al. Mathematical models for the aedes aegypti dispersal dynamics: travelling waves by wing and wind. Bulletin of Matematical Biology, 67 (2005), 509-528.
• ²⁵
  C. Tsallis. Introduction to Nonextensive statistical mechanics. Springer, (2009).
• ²⁶
  C. Tsallis & D.J. Bukman. Anomalous diffusion in the presence of external forces: Exact time-dependent solutions and their thermostatistical basis. Phys. Rev. E, 54 (1996), R2197.
• ²⁷
  S. Umarov, C. Tsallis & S. Steinberg. On a q-central limit theorem consistent with nonextensive statistical mechanics. Milan Journal of Mathematics, 76 (2008), 307-328.

• †
  Trabalho financiado pelo Cefet-MG e apresentado no CNMAC 2016.

Publication Dates

• Publication in this collection
  Aug 2017

History

• Received
  23 Nov 2016
• Accepted
  17 Apr 2017