In 1897, Biernacki introduced a blood test, known as Erythrocyte Sedimentation Rate (ESR), which helped in diagnosing the acute phase of inflammatory diseases and in following up the inflammatory process itself ^[2,3,4]. The discovery was announced in two articles ^[5,6]. At the beginning of nineteenth century, Fahraeus and Westergren, when performing pregnancy and tuberculosis tests, developed a test similar to ESR known as the Fahraeus-Westergren test ^[7,8,9,10].

Nowadays, due to the discovery of new and more accurate tests, ESR is little used despite its being a quick and low cost test. Nevertheless, the test is still recommended for patients with suspected giant cell arteritis, rheumatic polymyalgia and rheumatoid arthritis, among others ^[11]. However, as ESR is not very specific, it is often necessary to conduct further tests in order to confirm the results obtained by means of ESR, in order to avoid false-positive and false-negative results which are likely to occur in the presence of factors whose influence on blood properties would affect the test's results ^[12,13,14], such as age, anemia and pregnancy, resulting in increased ESR; polycythemia and increased leukocyte counting, resulting in decreased ESR; and analytic factors such as an inclined tube and room temperature, which would respectively increase and decrease ESR ^[15]. Other factors which affect the results are the presence of external vibration and tube deformation ^[16].

The concentration of nutrients in blood also plays a role in the analysis of ESR results ^[17]. Moreover, Nayha ^[18] noted that people who drink coffee and smoke present higher values of ESR. The use of some types of anticoagulants such as sodium citrate, oxalate or $K_{3}$ EDTA can also influence test results ^{[19,20,21,22]}.

Whelan et al. ^[23] published a work in which they measured the concentration of red cells at different times in blood samples of 5 male donors. In the same year, Huang et al. ^[24] developed a mathematical model to describe the behavior of the concentration of blood cells. Another notable work in ESR context was written by Sartory ^[25], whose aim was to study the prediction of erythrocyte sedimentation profiles. Moved by Huang's 1971 work, in 1990 Reuben and Shannon ^[26] discussed some problems in the mathematical modeling of the concentration of red blood cells. However, the authors of those studies did not take into account the transfer of nutrients from capillaries to tissues. Due to this fact, Sharma et al. ^[1] established a more precise mathematical model which takes into acoount such transfers.

The ESR test can be studied as a particular type of transport phenomenon ^[27]. It is worth mentioning that there exist several transport phenomena whose fractional models provide better descriptions than the corresponding classical models ^[28,29,30].

Our goal in studying the concentration of nutrients in blood is to show how fractional calculus employing a derivative in the Caputo sense provides a more realistic model in comparison to the classical one, i.e., the model with an integer order derivative.

In this work we assume an average speed equal to zero, thus restricting ourselves to the diffusion case. We use this model to introduce the basic concepts of fractional calculus and to present our fractional mathematical model. We propose a model with fractional derivatives in the Caputo sense with a time derivative of order $0 < \mu\leq 1$.

The solution obtained for the fractional mathematical model is given in terms of the Mittag-Leffler function and the Wright function. The solution has an extra degree of freedom in parameter $\mu$ ($0 < \mu \leq$ 1), which allows for a better fitting of experimental data on nutrient concentration in blood.

This paper is organized as follows: In section one we introduce the so-called fractional mathematical model associated with ESR, a generalization of the model proposed by Sharma et al. ^[1], which will be recovered through a limit process. Section two, our main result, is dedicated to obtaining the analytic solution of our model, which is found using the Laplace transform method and is expressed in terms of the Mittag-Leffler function and the Wright function. We also present a graphical analysis of the solution. In section three we recover as a special case, through an adequate limit process, the solution found by Sharma et al. ^[1]. Concluding remarks close the paper.

---

1. Time-fractional partial differential equation

The mathematical model proposed by Sharma et al. ^[1] describes the concentration of nutrients in blood by means of a non-homogeneous linear convection-diffusion partial differential equation (PDE). In this section we present a fractional version of that linear PDE. We assume that the average fluid velocity is equal to zero, i.e., we restrict our study to the diffusion case ^[31]. Our model can be considered a generalization of the Sharma et al. ^[1] model, in the sense that it recovers the latter as a special case, as we shall see in section three.

In this model, the concentration of nutrients in blood is a function $C(x, t)$ twice continuously differentiable that satisfies the following non-homogeneous time-fractional PDE,

$D_{L}\mathcal{D}_{x}^{2}C\left( x, t\right) - \mathcal{D}_{t}^{\mu }C\left( x, t\right) = \phi \left( x, t\right),$

(1)

with $0 < \mu\leq 1$, where $D_{L}$ is a positive constant and $\phi(x, t)$ is a twice continuously differentiable function describing the nutrient transfer rate and which satisfies the PDE

$D\mathcal{D}_{x}^{2}\phi \left( x, t\right) -k\phi \left( x, t\right) -% \mathcal{D}_{t}\phi \left( x, t\right) = 0,$

(2)

with both $D$ and $k$ positive constants.

The initial and boundary conditions imposed are given by

$\begin{equation*} \left\{ \begin{array}{ll} \phi (x, 0) = \exp \left( -\sqrt{\frac{k-a}{D}}x\right), & k\geq a, D>0, \\ \phi (0, t) = \exp \left( -at\right), & t>0, \\ \phi (\infty, t) = 0, & t>0. \end{array}% \right. \end{equation*}$

The solutions of Eq.(2) can be written as

$\phi \left( x, t\right) = \exp \left( -\left( at+bx\right) \right),$

where $b^{2} = \displaystyle\frac{\left(k-a\right) }{D}>0$ and $a$ is a constant to be adequately chosen from a known value of $\phi \left(x, t\right)$.

We assume that the fractional derivative of order $\mu$, $0 < \mu \leq 1$ is considered in the Caputo sense ^[32,33,34], defined as follows:

$\begin{equation*} \mathcal{D}_{t}^{\mu }C\left( x, t\right) : = \left\{ \begin{array}{c} \displaystyle\frac{1}{\Gamma \left( n-\mu \right) }\int_{0}^{t}C^{\left( n\right) }\left( \tau, t\right) \left( t-\tau \right) ^{n-\mu -1}d\tau, \;\;n-1<\mu <n \\ \\ C^{\left( n\right) }\left( x, t\right), \quad \quad \quad \quad \quad \quad \quad \quad \mu = n, \end{array} \right. \end{equation*}$

where $\mathcal{D}^{\mu}_t \equiv \displaystyle\frac{\partial^{\mu }}{\partial t^{\mu}}$ and $C^{(n)}(x, t)$ is the usual derivative of order $n$ with respect to $t$, $C^{n}(x, t)\in AC^{n}[0, h]$, where $AC^{n}[0, h]$ is the space of absolutely continuous functions and $t>0$. Furthermore, we must impose the following initial and boundary conditions for Eq.(1):

$\begin{equation} \left\{ \begin{array}{ll} C(x, 0) = 0, & x\geq 0 \\ C(0, t) = 1, & t>0 \\ C(\infty, t) = 0, & t>0, \end{array}% \right. \end{equation}$

(3)

with $C(x, t)\in C^{2}[0, h]$.

$AC^{n}\left[a, b\right] = \left\{ f:\left[a, b\right] \rightarrow \mathbb{C} \text{ and }\left(D^{n-1}f\right) \left(x\right) \in AC\left[a, b\right] \text{ where }\left(D = \frac{d}{dx}\right) \right\}$.

Thus, from these considerations, it follows that the time-fractional mathematical model to be addressed is composed of a non-homogeneous fractional PDE

$D_{L}\mathcal{D}_{x}^{2}C\left( x, t\right) - \mathcal{D}_{t}^{\mu }C\left( x, t\right) = \exp \left( -\left( at+bx\right) \right), \quad a, b \in \mathbb{R},$

(4)

with initial and boundary conditions given by Eq.(3).

---

2. Analytic solution

In this section, we solve this problem, employing the methodology of Laplace transform to convert the non-homogeneous fractional PDE into a non-homogeneous linear ordinary differential equation.

Then, applying the Laplace transform ^[35,36] in the time variable $t$ on both sides of Eq.(4), we have

$D_{L}\frac{d^{2}}{dx^{2}}C\left( x, s\right) -s^{\mu }C\left( x, s\right) +s^{\mu -1}C\left( x, 0\right) = \frac{\exp \left( -bx\right) }{s+a}.$

Using the initial condition $C(x, 0) = 0$ we can rewrite this equation as

$D_{L}\frac{d^{2}}{dx^{2}}C\left( x, s\right) -s^{\mu }C\left( x, s\right) =\frac{\exp \left( -bx\right) }{s+a},$

(5)

where $0 < \mu \leq 1$, $D_{L}>0$ and

$C\left( x, s\right) =\mathscr{L}\left\{ C\left( x, t\right) \right\} = \colon \int_{0}^{\infty }e^{-st}C\left( x, t\right) dt$

is the Laplace transform of $C(x, t)$ with parameter $s$, $\mbox{Re}(s)>0$. We assume that $C(x, t)$ is continuous by parts on $[0, \infty]$ and of exponential order.

Using the methods of characteristic equation and undetermined coefficients in Eq.(5) we obtain the general solution, given by

$C\left( x, s\right) = \left( \frac{1}{s}+\frac{1}{\left( s+a\right) \left( s^{\mu }-\frac{b^{2}}{\alpha ^{2}}\right) }\right) \exp \left( -\alpha xs^{\mu /2}\right) +\frac{\exp \left( -bx\right) }{\left( s+a\right) \left( \frac{b^{2}}{\alpha ^{2}}-s^{\mu }\right) },$

(6)

where $\alpha ^{2} = \displaystyle\frac{1}{D_{L}}$ and $D_{L}>0$.

In order to recover the solution in the time variable we take the inverse Laplace transform on both sides of Eq.(6), obtaining

$\begin{eqnarray} \label{S16} C\left( x, t\right) & = &\mathscr{L}^{-1}\left\{ C\left( x, s\right) \right\} = % \mathscr{L}^{-1}\left\{ \frac{1}{s}\exp \left( -\alpha xs^{\mu /2}\right) \right\} + \notag \\ &&+\mathscr{L}^{-1}\left\{ \frac{1}{\left( s+a\right) \left( s^{\mu }-\frac{% b^{2}}{\alpha ^{2}}\right) }\exp \left( -\alpha xs^{\mu /2}\right) \right\} \notag \\ &&-\mathscr{L}^{-1}\left\{ \frac{1}{\left( s+a\right) \left( s^{\mu }-\frac{% b^{2}}{\alpha ^{2}}\right) }\exp \left( -bx\right) \right\}, \end{eqnarray}$

(7)

where

$C\left( x, t\right) =\mathscr{L}^{-1}\left\{ C\left( x, s\right) \right\} = \colon \frac{1% }{2\pi i}\int_{\gamma -i\infty }^{\gamma +i\infty }e^{st}C\left( x, s\right) ds$

is the inverse Laplace transform and the integral is performed in the complex plane with the singularities $C(x, s)$ on the left side of $\gamma =\mbox{Re}(s)$ ^[31].

Introducing the change $\beta ^{2} = b^{2}D_{L}$, we rewrite Eq.(7) as

$C\left( x, t\right) = C_{1}\left( x, t\right) +C_{2}\left( x, t\right) - \mbox{exp}(-bx) C_{3}\left( x, t\right),$

with

$C_{1}\left(x, t\right) = \mathscr{L}^{-1}\left\{ \frac{\exp \left( -\alpha xs^{\mu /2}\right) }{s}\right\} ;$

(8)

$C_{2}\left( x, t\right) = \mathscr{L}^{-1}\left\{ \frac{\exp \left( -\alpha xs^{\mu /2}\right) }{\left( s+a\right) \left( s^{\mu }-\beta ^{2}\right) }% \right\} ;$

(9)

$C_{3}\left( x, t\right) = \underset{x\rightarrow 0}{\lim }C_{2}\left( x, t\right) .$

We then calculate each inverse Laplace transform separately. To calculate $C_{1}\left(x, t\right)$ we introduce the MacLaurin series associated with the exponential function; choosing $f^{(k)}(0) = 1$ in the series, we have

$\frac{1}{s}\exp \left( -\alpha xs^{\mu /2}\right) = \sum\limits_{k = 0}^\infty\frac{\left( -\alpha x\right) ^{k}}{k!}s^{\frac{\mu k}{ 2}-1}.$

(10)

Applying the inverse Laplace transform on both sides of Eq.(10), and using the result

$\mathscr{L}^{-1}\left\{ s^{-q}\right\} = \frac{t^{q-1}}{\Gamma \left( q\right)},$

with $\mbox{Re}(q)>0$, $q = 1-\mu k/2$, we can rewrite Eq.(8) as follows:

$C_{1}\left( x, t\right) = \sum\limits_{k = 0}^\infty \frac{\left( -\alpha x/t^{\mu /2}\right) ^{k}}{k!\Gamma \left( 1-\mu k/2\right) }.$

(11)

Moreover, considering $\beta = 1$, $\alpha = -\mu/2$ and $z = -\displaystyle\frac{\alpha x}{t^{\mu /2}}$ we obtain

$C_{1}\left( x, t\right) = \mathbb{W}\left( -\mu /2, 1;-\frac{\alpha x}{t^{\mu /2}}\right) .$

(12)

where

$\mathbb{W}\left( -\mu /2, 1;z\right) = \sum\limits_{k = 0}^\infty \frac{z^{k}}{k!\Gamma \left( -\mu k/2+1\right) }.$

(13)

is the Wright function ^[37].

We now evaluate the second inverse Laplace transform. As with $C_{1}(x, t)$, we also write the exponential function in terms of its MacLaurin series. Once more, applying the inverse Laplace transform we can write

$\mathscr{L}^{-1}\left\{ \frac{1}{\left( s+a\right) \left( s^{\mu }-\beta ^{2}\right) } \exp \left( -\alpha xs^{\mu /2}\right) \right\} = \sum\limits_{m = 0}^\infty\frac{\left( -\alpha x\right) ^{m}}{m!}% \mathscr{L}^{-1}\left\{ \frac{s^{\mu m/2}}{\left( s+a\right) \left( s^{\mu }-\beta ^{2}\right) }\right\} .$

(14)

In order to evaluate this inverse Laplace transform, we consider the following expression ^[38]:

$\Omega = \displaystyle\frac{s^{\sigma }}{s^{\alpha }+\widetilde{a}s^{\delta }+bs^{\gamma }+cs^{\mu }+d},$

with $\widetilde{a}, b, c, d\in\mathbb{R}$ and $\alpha, \delta, \gamma, \mu\in\mathbb{R}$ such that $\widetilde{a}\neq 0$ and $\alpha>\delta>\gamma>\mu$.

Assuming the condition $\left\vert \displaystyle\frac{bs^{\gamma }+cs^{\mu }+d}{s^{\alpha }+\widetilde{a}s^{\delta }}\right\vert < 1$ and using the geometric series we have

$\begin{eqnarray}\label{S28} \sum\limits_{k = 0}^\infty \left( -1\right) ^{k}s^{\sigma -\delta-\delta k}\frac{\left( bs^{\gamma }+cs^{\mu }+d\right) ^{k}}{\left( s^{\alpha -\delta }+\widetilde{a}\right) ^{k+1}} & = &\frac{s^{\sigma }}{s^{\alpha }+s^{\delta }\widetilde{a}}\left( \frac{1}{1+\displaystyle\frac{bs^{\gamma }+cs^{\mu }+d}{s^{\alpha }+\widetilde{a}s^{\delta }}}\right) \notag \\ & = &\frac{s^{\sigma }}{bs^{\gamma }+cs^{\mu }+d+s^{\alpha }+\widetilde{a}s^{\delta }}. \end{eqnarray}$

(15)

The binomial theorem and the definition of binomial coefficients ^[39] allow us to rewrite Eq.(15) as

$\begin{eqnarray}\label{S29} \Omega & = &\sum\limits_{k = 0}^\infty \left( -1\right) ^{k}% \sum\limits_{l = 0}^k \binom{k}{l}d^{l}\left( bs^{\gamma }+cs^{\mu }\right) ^{k-l}\frac{s^{\sigma -\delta -\delta k}}{\left( s^{\alpha -\delta }+\widetilde{a}\right) ^{k+1}} \notag \\ & = &\sum\limits_{k = 0}^\infty \left( -1\right) ^{k}\overset{k}{% \underset{l = 0}{\sum }}\frac{k!}{l!\left( k-l\right) !}d^{l}\underset{j = 0}{% \overset{k-l}{\sum }}\frac{\left( k-l\right) !}{j!\left( k-l-j\right) !}% \left( bs^{\gamma }\right) ^{k-l-j}\left( cs^{\mu }\right) ^{j}\frac{% s^{\sigma -\delta -\delta k}}{\left( s^{\alpha -\delta }+\widetilde{a}\right) ^{k+1}} \notag \\ & = &\sum\limits_{k = 0}^\infty \left( -1\right) ^{k}b^{k}k!% \sum\limits_{l = 0}^k \frac{\left( d/b\right) ^{l}}{l!}\underset{ j = 0}{\overset{k-l}{\sum }}\frac{\left( c/b\right) ^{j}}{j!\left( k-l-j\right) !}\Lambda _{\sigma }, \end{eqnarray}$

(16)

where $\Lambda _{\sigma } = \displaystyle\frac{s^{\sigma -\delta \left(1+k\right) +\mu j+\gamma\left(k-l-j\right) }}{\left(s^{\alpha -\delta }+\widetilde{a}\right) ^{k+1}}$.

Taking the inverse Laplace transform on both sides of Eq.(16) and using the result

$\mathscr{L}^{-1}\left\{ \Lambda _{\sigma }\right\} = \mathscr{L}^{-1}\left\{ \frac{s^{\sigma -\delta \left( 1+k\right) +\mu j+\gamma \left( k-l-j\right) }% }{\left( s^{\alpha -\delta }+\widetilde{a}\right) ^{k+1}}\right\} = t^{\xi -1}\mathbb{E}% _{\alpha -\delta, \xi }^{k+1}\left( -\widetilde{a}t^{\alpha -\delta }\right),$

(17)

we get

$\mathscr{L}^{-1}\left\{ \Omega \right\} =\sum\limits_{k = 0}^\infty\left( -1\right) ^{k}b^{k}k!\sum\limits_{l = 0}^k \frac{% \left(d/b\right) ^{l}}{l!}\sum\limits_{j = 0}^{k-l} \frac{% \left(c/b\right) ^{j}}{j!\left( k-l-j\right) !}t^{\xi -1} \mathbb{E}_{\alpha -\delta, \xi}^{k+1}\left( -\widetilde{a}t^{\alpha -\delta }\right),$

(18)

with $\xi = -\sigma +\alpha +\left(\alpha -\gamma \right) k+\gamma l-\left(\mu-\gamma \right) j$ and where $\mathbb{E}_{\alpha -\delta, \xi}^{k+1}\left(\cdot\right)$ is the three-parameters Mittag-Leffler function ^[38,40].

In particular, considering $c = 0$ in Eq.(18), we have that $j = 0$ is the only term contributing to the sum and we conclude that

$\mathscr{L}^{-1}\left\{ \frac{s^{\sigma }}{s^{\alpha }+\widetilde{a}s^{\delta }+bs^{\gamma }+d}\right\} = \sum\limits_{k = 0}^\infty \left( -1\right) ^{k}b^{k}k!\sum\limits_{l = 0}^k \frac{\left( d/b\right) ^{l}t^{\xi -1}}{l!\left( k-l\right) !}\mathbb{E}_{\alpha -\delta , \xi }^{k+1}\left( -\widetilde{a}t^{\alpha -\delta }\right),$

(19)

where $\xi = -\sigma +\alpha +\left(\alpha -\gamma \right) k+\gamma l$ and $\alpha >\delta >\gamma$.

Then, putting $\sigma = \mu m/2$, $d = -a\beta ^{2}$, $\alpha = \mu +1$, $\gamma = \mu$, $\delta = 1$, $b = a$ and $\widetilde{a} = -\beta ^{2}$ in Eq.(19) and going back to Eq.(14), we can write

$C_{2}\left( x, t\right) = t^{\mu }\sum\limits_{m = 0}^\infty% \frac{\left( -\alpha xt^{-\mu /2}\right) ^{m}}{m!}\sum\limits_{k = 0}^\infty\left( -at\right) ^{k}k!\sum\limits_{l = 0}^k % \frac{\left( -\beta ^{2}t^{\mu }\right) ^{l}}{l!\left( k-l\right) !}\mathbb{E% }_{\mu, \theta }^{k+1}\left( \beta ^{2}t^{\mu }\right),$

(20)

where $\theta = -\mu m/2+\mu +1+k+\mu l$.

In order to write the solution of the PDE in terms of the two-parameters Mittag-Leffler function, we evaluated the sum on $l$ appearing in the last expression in order to find a relationship between two-and three-parameters Mittag-Leffler functions. Using the identity

$\Lambda = \sum\limits_{j = 0}^k\frac{\left( z\right) ^{j}}{% j!\left( k-j\right) !}\mathbb{E}_{\lambda, \lambda j+\delta }^{\rho }\left( -z\right) = \sum\limits_{j = 0}^k\sum\limits_{l = 0}^\infty \frac{\left( z\right) ^{j}}{j!\left( k-j\right) !}\frac{\left( \rho \right)_{l}\left( -z\right) ^{l}}{l!\Gamma \left( \lambda l+\lambda j+\delta \right) },$

(21)

where $\left(\rho \right) _{l} = \rho \left(\rho +1\right)....\left(\rho+l-1\right)$, together with the definition and properties of the binomial coefficients in Eq.(21), we can write ^[38]

$\begin{eqnarray}\label{S40} \sum\limits_{j = 0}^k\frac{\left( z\right) ^{j}}{j!\left( k-j\right) !}\mathbb{E}_{\lambda, \lambda j+\delta }^{\rho }\left( -z\right) & = &\underset{i = 0}{\overset{\infty }{\sum }}\frac{\left( -z\right) ^{i}}{% \Gamma \left( \lambda i+\delta \right) }\frac{1}{k!}\sum\limits_{j = 0}^k \frac{\left( -1\right) ^{j}k!}{j!\left( k-j\right) !}\binom{i-j+\rho -1}{\rho -1} \notag \\ & = &\underset{i = 0}{\overset{\infty }{\sum }}\frac{\left( -z\right) ^{i}}{% \Gamma \left( \lambda i+\delta \right) }\frac{1}{k!}\frac{\left( \rho -k\right) _{i}}{i!} = \frac{1}{k!}\mathbb{E}_{\lambda, \delta }^{\rho -k}\left( -z\right) . \end{eqnarray}$

(22)

Choosing $z = -\beta ^{2}t^{\mu }$, $\rho = k+1$, $\lambda = \mu$, $j = l$ and $\delta = k+\mu +1-\mu m/2$ in Eq.(22) and substituting the result into Eq.(20), we conclude that

$C_{2}\left( x, t\right) = t^{\mu }\sum\limits_{m = 0}^\infty% \frac{\left( -\alpha xt^{-\mu /2}\right)^{m} }{m!}\sum\limits_{k = 0}^\infty \left( -at\right) ^{k}\mathbb{E}_{\mu, \mu +k+1-\mu m/2}\left( \beta ^{2}t^{\mu }\right),$

(23)

where $\mathbb{E}_{\alpha, \beta }\left(\cdot\right)$ is the two-parameters Mittag-Leffler function, which is considered uniformely convergent ^[40].

The last inverse Laplace transform, $C_{3}(x, t)$, is obtained by means of an adequate limit, i.e., we consider $x\rightarrow 0$ in Eq.(23). The only term that contributes in this limit is $m = 0$, i.e., we get

$C_{3}\left( x, t\right) = t^{\mu }\sum\limits_{k = 0}^\infty % \left( -at\right) ^{k}\mathbb{E}_{\mu, \mu +k+1}\left( \beta ^{2}t^{\mu }\right) .$

(24)

Thus, from the results obtained in Eq.(12), Eq.(23) and Eq.(24), we get the solution of our initial problem, i.e., a solution of Eq.(4) satisfying the conditions given by Eq.(3):

$\begin{eqnarray} \label{S43} C\left( x, t\right) & = &t^{\mu }\sum\limits_{m = 0}^\infty\frac{% \left( -\alpha xt^{-\mu /2}\right) ^{m}}{m!}\underset{k = 0}{\overset{\infty }{% \sum }}\left( -at\right) ^{k}\mathbb{E}_{\mu, \mu +k+1-\mu m/2}\left( \beta ^{2}t^{\mu }\right) + \\ &&+\mathbb{W}\left( -\mu /2, 1;-\frac{\alpha x}{t^{\mu /2}}\right) -\exp \left( -bx\right) t^{\mu }\sum\limits_{k = 0}^\infty \left( -at\right) ^{k}\mathbb{E}_{\mu, \mu +k+1}\left( \beta ^{2}t^{\mu }\right) \notag, \end{eqnarray}$

(25)

where the parameters are given by $\alpha ^{2} = 1/D_{L}$, $\beta ^{2} = b^{2}D_{L}$ and $0 < \mu \leq 1$. The solution given by Eq.(25}) valid for $t>0$ is $AC^{n}[0, h]$ and class $C^{2}[0, h]$; then, substituting it into Eq.(4) we can easily verify that it satisfies the IVP (Initial value problem) and BVP (Boundary value problem) [Eq.(2) and Eq.(3)] ^[31].

Let us now perform a graphical analysis. For this sake, we have to choose values for some parameters appearing in the solution given by Eq.(25). We used the following values: axial dispersion coefficient $D_{L} =4.8\times 10^{-4}\mbox{cm}^2 \mbox{s}^{-1}$ ^[41]; diffusivity coefficient of oxygen $D = 9.8\times 10^{-5}\mbox{cm}^2 \mbox{s}^{-1}$ ^[42]; nutrient transfer coefficient $k = 1.5\times 10^{-4}\mbox{m}\mbox{s}^{-1}$ ^[43]; $a = -0.005\times 10^{-4}\mbox{m}\mbox{s}^{-1}$ ^[43]. We also fix a time $t = 15\mbox{s}$ and we consider the interval $x =[0, 4]$ (which can be extended).

In Figures 1 and 2, the horizontal axis $x$ represents space and the vertical axis $y$ is the normalized concentration of nutrients in blood.

The parameter values used to plot Figure 1 were also used to plot the solution of the integer order PDE, Figure 2. The graphics were plotted using MATLAB 7:10 software (R2010a).

Remark that as $x$ (space) increases, the value of $C/C_{1}$ (concentration of nutrients) decreases, that is, when we move towards the extremity of the artery ($x\neq 0$), the blood concentration of solute decreases. A decrease in solute concentration means that cells are not enough efficient in getting their nutrition, so we conclude that the efficiency of nutrient transport near the artery is greater than at its venous extremity.

As we have already said, with the freedom provided by parameter $\mu$ ($0 < \mu\leq1$), it is possible to describe more accurately the information about the concentration of nutrients near the arterial extremity because, as seen above, the fractionalization of the derivative refines the solution. Note that for $\mu =0.10$ the behavior of the analytic solution remains near the arterial ($x =0$) for longer time. We can thus see that as $\mu\rightarrow 1$, the fractional solution converges to the solution of the integer order PDE.

We supposed that the space variable $x$ lies within the range $[0, 4]$. We might as well have analyzed variable $x$ in the range $[0, 12]$ or any other interval; however, the first representative interval is the one we chose because for $x\geq3.8$ the level $C/C_{1}$ remains below the $x$ axis. So it is interesting, in this context, to carry our analysis only on the $[0, 4]$ range.

---

3. Particular case: $\mu \rightarrow 1$

In this section, we analyze the solution of the fractional PDE in the limit $\mu \rightarrow 1$, in order to recover the result found by Sharma et al. ^[1].

Since the solution of the fractional PDE Eq.(4) is given by Eq.(25), taking the limit $\mu\rightarrow 1$, it follows that

$\begin{eqnarray} \label{S44} C\left( x, t\right) & = &t\sum\limits_{m = 0}^\infty\frac{\left( -\alpha xt^{-1/2}\right) ^{m}}{m!}\sum\limits_{k = 0}^\infty % \left( -at\right) ^{k}\mathbb{E}_{1, k+2-m/2}\left( \beta ^{2}t\right)+ \notag \\ &&+\mathbb{W}\left( -1/2, 1;-\frac{\alpha x}{t^{1/2}}\right) -\exp \left( -bx\right) t\sum\limits_{k = 0}^\infty \left( -at\right)^{k}% \mathbb{E}_{1, k+2}\left( \beta ^{2}t\right) . \end{eqnarray}$

(26)

In the last two terms of the sum in Eq.(26), we can use the results involving the Wright function and the complementary error function and the exponential function, to get ^[37]:

$\begin{eqnarray} \label{S47} C\left( x, t\right) & = &t\sum\limits_{m = 0}^\infty\frac{\left( -\alpha xt^{-1/2}\right) ^{m}}{m!}\sum\limits_{k = 0}^\infty % \left( -at\right) ^{k}\mathbb{E}_{1, k+2-m/2}\left( \beta ^{2}t\right) + \notag \\ &&+1+\mbox{erf}\left( -\alpha x/2\sqrt{t}\right) -\exp \left( -bx\right) \frac{\exp \left( \beta ^{2}t\right) -\exp \left( -at\right) }{a+\beta ^{2}}. \end{eqnarray}$

(27)

We want to express Eq.(27) in terms of $\mbox{erfc}(\cdot)$ and $\mbox{exp}(\cdot)$. We then evaluate the inverse Laplace transform in Eq.(9) using partial fractions.

Taking the limit $\mu \rightarrow 1$ in Eq.(9), it follows that

$C_{2}\left( x, t\right) = \mathscr{L}^{-1}\left\{ \frac{\exp \left( -\alpha x% \sqrt{s}\right) }{\left( s+a\right) \left( s-\beta ^{2}\right) }\right\} .$

(28)

Using partial fractions and taking the inverse Laplace transform, we have

$\begin{eqnarray} \label{S51} 2\left( \beta ^{2}+a\right) \mathscr{L}^{-1}\left\{ \frac{\exp \left( -\alpha x\sqrt{s}\right) }{\left( s+a\right) \left( s-\beta ^{2}\right) }% \right\} & = &-\mathscr{L}^{-1}\left\{ \frac{\exp \left( -\alpha x\sqrt{s}% \right) }{\sqrt{s}\left( \sqrt{s}-i\sqrt{a}\right) }\right\} -\mathscr{L}% ^{-1}\left\{ \frac{\exp \left( -\alpha x\sqrt{s}\right) }{\sqrt{s}\left( \sqrt{s}+i\sqrt{a}\right) }\right\} + \notag \\ &&+\mathscr{L}^{-1}\left\{ \frac{\exp \left( -\alpha x\sqrt{s}\right) }{% \sqrt{s}\left( \sqrt{s}-\beta \right) }\right\} +\mathscr{L}^{-1}\left\{ \frac{\exp \left( -\alpha x\sqrt{s}\right) }{\sqrt{s}\left( \sqrt{s}+\beta \right) }\right\} . \end{eqnarray}$

(29)

In evaluating the inverse Laplace transforms, we can use the following result ^[44]:

$\mathscr{L}^{-1}\left\{ \frac{\exp \left( -k\sqrt{s}\right) }{\sqrt{s}\left( \sqrt{s}+b\right) }\right\} = \exp \left( bk\right) \exp \left( b^{2}t\right) % \mbox{erfc}\left( b\sqrt{t}+\frac{k}{2\sqrt{t}}\right),$

(30)

with $k\geq 0$, $b\in \mathbb{C}$ and where $\mbox{erfc}\left(x\right)$ is the complementary error function.

Thus, applying Eq.(30) in each term of Eq.(29), we have

$\begin{eqnarray}\label{S54} 2\left( \beta ^{2}+a\right) \mathscr{L}^{-1}\left\{ \frac{1}{\left( s+a\right) \left( s-\beta ^{2}\right) }\right\} & = &\exp \left( \beta ^{2}t\right) \left( \mbox{erfc}\left( \beta \sqrt{t}\right) +\mbox{erfc}\left( -\beta \sqrt{t}\right) \right) - \notag \\ &&-\exp \left( -at\right) \left( \mbox{erfc}\left( i\sqrt{at}\right) +% \mbox{erfc}\left( -i\sqrt{at}\right) \right).\notag\\ \end{eqnarray}$

(31)

Analyzing the error functions in Eq.(31), we conclude that

$\mathscr{L}^{-1}\left\{ \frac{1}{\left( s+a\right) \left( s-\beta ^{2}\right) }\right\} = \frac{\exp \left( \beta ^{2}t\right) -\exp \left( -at\right) }{\beta ^{2}+a}.$

(32)

As we evaluated the inverse Laplace transform of $C_ {2}(x, t)$ in Eq.(9) in the case $\mu =1$ using two different procedures, involving respectively a Mittag-Leffler function and error functions, we can write, as a by-product, the following interesting mathematical identity involving Mittag-Leffler functions:

$\begin{eqnarray} \label{S56} &&2\left( a+\beta ^{2}\right) t\sum\limits_{m = 0}^\infty\frac{% \left( -\alpha xt^{-1/2}\right) }{m!}\sum\limits_{k = 0}^\infty % \left( -at\right) ^{k}\mathbb{E}_{1, 2+k-m/2}\left( \beta ^{2}t\right) \begin{array}{c} =\end{array} \notag \\ & = &e^{\beta \alpha x}e^{\beta ^{2}t}\mbox{erfc}\left( \beta \sqrt{t}+\frac{% \alpha x}{2\sqrt{t}}\right) +e^{-\beta \alpha x}e^{\beta ^{2}t}\mbox{erfc}% \left( -\beta \sqrt{t}+\frac{\alpha x}{2\sqrt{t}}\right) - \notag \\ &&-e^{i\alpha \sqrt{a}x}e^{-at}\mbox{erfc}\left( i\sqrt{at}+\frac{\alpha x}{2% \sqrt{t}}\right) -e^{-i\alpha \sqrt{a}x}e^{-at}\mbox{erfc}\left( -i\sqrt{at}+% \frac{\alpha x}{2\sqrt{t}}\right). \notag\\ \end{eqnarray}$

(33)

Further, considering $\alpha = 0$ in Eq.(33), which means that only $m = 0$ contributes to the first sum, we obtain

$t\sum\limits_{k = 0}^\infty \left( -at\right) ^{k}\mathbb{E}% _{1, 2+k}\left( \beta ^{2}t\right) = \frac{e^{\beta ^{2}t}-e^{-at}}{\left( a+\beta ^{2}\right) }.$

(34)

Consequently, Eq.(33) can be interpreted as a generalization of Eq.(34). Also, considering $a = 0$ in the previous equation, we have

$\beta ^{2}t\mathbb{E}_{1, 2}\left( \beta ^{2}t\right) = e^{\beta ^{2}t}-1,$

which is a known identity involving the Mittag-Leffler function ^[40].

Finally, we can write the main relation we need to recover the solution proposed by Sharma et al. ^[1]. According to Eq.(27) and Eq.(33):

$\begin{eqnarray}\label{S57} C\left( x, t\right) & = &t\sum\limits_{m = 0}^\infty\frac{\left( -\alpha xt^{-1/2}\right) ^{m}}{m!}\sum\limits_{k = 0}^\infty % \left( -at\right) ^{k}\mathbb{E}_{1, k+2-m/2}\left( \beta ^{2}t\right) + \notag \\ &&+1+\mbox{erf}\left( -\alpha x/2\sqrt{t}\right) -\exp \left( -bx\right) \frac{\exp \left( \beta ^{2}t\right) -\exp \left( -at\right) }{a+\beta ^{2}} \notag \\ & = &1-\mbox{erf}\left( \alpha x/2\sqrt{t}\right) -\frac{\exp \left( -bx\right) }{a+\beta ^{2}}\left( \exp \left( \beta ^{2}t\right) -\exp \left( -at\right) \right) \notag \\ &&+\frac{\exp \left( \beta ^{2}t\right) }{2\left( a+\beta ^{2}\right) }\left[ \begin{array}{c} \exp \left( \beta \alpha x\right) \mbox{erfc}\left( \beta \sqrt{t}+\frac{% \alpha x}{2\sqrt{t}}\right) \\ +\exp \left(-\beta \alpha x\right) \mbox{erfc}\left(-\beta \sqrt{t}+\frac{% \alpha x}{2\sqrt{t}}\right) \end{array}% \right] - \notag \\ &&-\frac{\exp \left( -at\right) }{2\left( a+\beta ^{2}\right) }\left[ \begin{array}{c} \exp \left( i\alpha \sqrt{a}x\right) \mbox{erfc}\left( i\sqrt{at}+\frac{% \alpha x}{2\sqrt{t}}\right) \\ +\exp \left(-i\alpha \sqrt{a}x\right) \mbox{erfc}\left(-i\sqrt{at}+\frac{% \alpha x}{2\sqrt{t}}\right) \end{array}% \right] . \end{eqnarray}$

(35)

We emphasize that parameters $D$ and $D_{L}$ are positive constants and $k\geq a$, as imposed in both models. Moreover, returning to the original parameters $\beta = \sqrt{\frac{\left(k-a\right) }{D}D_{L}}$, $b = \sqrt{\frac{k-a}{D}}$, from Eq.(35), we conclude that

$\begin{eqnarray}\label{S58} C\left( x, t\right) & = &\mbox{erf}c\left( \frac{x}{2\sqrt{D_{L}t}}\right) -% \frac{\exp \left( -\sqrt{\frac{k-a}{D}}x\right) }{k\left( \frac{D_{L}}{D}% \right) +a\left( 1-\frac{D_{L}}{D}\right) }\left( \exp \left( \left( \frac{% k-a}{D}\right) D_{L}t\right) -\exp \left( -at\right) \right) + \notag \\ &&+\frac{\exp \left( \left( \frac{k-a}{D}\right) D_{L}t\right) }{2\left[ k\left( \frac{D_{L}}{D}\right) +a\left( 1-\frac{D_{L}}{D}\right) \right] }% \left[ \begin{array}{c} \exp \left( \sqrt{\frac{k-a}{D}}x\right) \mbox{erfc}\left( \frac{x+2D_{L}t% \sqrt{\frac{k-a}{D}}}{2\sqrt{D_{L}t}}\right) \\ +\exp \left(-\sqrt{\frac{k-a}{D}}x\right) \mbox{erfc}\left( \frac{x-2D_{L}t% \sqrt{\frac{k-a}{D}}}{2\sqrt{D_{L}t}}\right) \end{array}% \right] - \notag \\ &&-\frac{\exp \left( -at\right) }{2\left[k\left( \frac{D_{L}}{D}\right) +a\left( 1-\frac{D_{L}}{D}\right) \right] }\left[ \begin{array}{c} \displaystyle\exp \left( \frac{i\sqrt{a}x}{\sqrt{D_{L}}}\right) \mbox{erfc}\left( \frac{% x+2it\sqrt{D_{L}a}}{2\sqrt{D_{L}t}}\right) \\ +\exp \left(-\frac{i\sqrt{a}x}{\sqrt{D_{L}}}\right) \mbox{erfc}\left( \frac{% x-2it\sqrt{D_{L}a}}{2\sqrt{D_{L}t}}\right) \end{array}% \right], \end{eqnarray}$

(36)

which is exactly the result obtained in ^[1].

---

4. Concluding remarks

After a brief introduction to the study of the concentration of nutrients in blood, a factor that interferes with ESR, we proposed a fractional mathematical model employing fractional derivatives in the Caputo sense. We obtained its analytic solution in terms of the Mittag-Leffler function and the Wright function using the methodology of Laplace transform in the time variable. We should point out that one of the greatest challenges of fractional calculus, in the study of differential equations, is to propose a fractional differential equation whose corresponding analytic solution recovers the integer order case in an adequate limit. Here, it was possible to recover the solution of the integer case applying the limit $\mu\rightarrow 1$ to the analytic solution, Eq.(25), of the fractional PDE, Eq.(4). As for what was expected about the relation between the fractional mathematical model and the integer order model of ^[1], we can say that our fractional model provides more accurate information about the concentration of nutrients in blood, as one can also see in Figure 1.

A natural continuation of this work is to confront our fractional model with experimental data, in order to be able to make predictions using ESR tests. Studies in this direction are being done and will be published in a forthcoming paper.

---

Acknowledgment

We thank Prof. Dr. Felix Silva Costa and Dr. J. Emílio Maiorino for fruitful discussions and Prof. Dr. F. Mainardi for suggesting several references on the subject. We also thank anonymous referees by the suggestions that improved the paper.

---

Conflict of Interest

All authors declare no conflicts of interest in this paper.

---