A Stochastic Approach to Modeling Food Pattern

Abstract

In this paper, we propose a fractional differential equation of order one-half, to model the evolution through time of the dynamics of accumulation and elimination of the contaminant in the human organism with a deficient immune system, during consecutive intakes of contaminated food. This process quantifies the exposure to toxins of subjects living with comorbidity (children not breastfed, the elderly, and pregnant women) to food-borne diseases. The Adomian Decomposition Method and the fractional integration of Riemann Liouville are used in the modeling processes.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Food safety is an important issue in the search for solutions to the various pathologies such as food-borne diseases, cancer, renal failure, diabetes, AIDS, Crohn’s disease, ulcerative colitis, and so on that the human species has been facing in recent decades. The continuous time model named the Kinetic Dietary Exposure Model (KDEM) introduced in [1] and later in [2] is a stochastic process which allows modeling the phenomena of accumulation and elimination of the contaminant from the human organism and which gives an answer to this quest for a solution to the current health problems. Their model aims to represent the evolution through time of a number of contaminants in the human body.

Note that this theory developed to model the evolution of contaminant dynamics in the body and does not take into account a failing immune system. Subjects who have contracted bacteria, malnourished children, diabetics etc., have a weak immune system. So, the speed of elimination of the contaminant in the body becomes slow. Note that the linear pharmacokinetic model with one compartment is used in the existing model to describe the phenomenon of elimination of the contaminant by the immune system. The quantity $Xt$ of the contaminant in the body between two consecutive intakes decreases according to the linear differential equation of the form:$$\begin{array}{}\text{(1)}& \frac{\mathrm{d}xt}{\mathrm{d}t}=-\theta xt,\end{array}$$where $Xt$ is given time t, $\theta >0$ is the removal rate of the contaminant between the dates $T_i$ and $T_{i+1},$ with $\theta =ln2/DV,$ DV being the biological half-life of the contaminant (see more in [3]).

The dynamics of the evolution of the accumulation and elimination of the contaminant proposed in the KDEM model is of the form:$$\begin{array}{}\text{(2)}& X_{n+1}t=X_n.e^{-\theta t}+E_{n+1},\end{array}$$and the exposure E of an individual is calculated for all p = 1, …, P as follows:$$\begin{array}{}\text{(3)}& E=\frac{{\displaystyle {\sum }_{p=1}^P{Q}_p{C}_p}}{w},\end{array}$$where P is the product, $Q_p$ is the contamination of p foods, $C_p$ is the consumption of p foods, and $w$ is the individual’s body weight (see[3]).

The structural properties of this exposure process have been studied extensively in [5] using integrated Markov chain analysis. Such autoregressive models with random coefficients have been widely studied in the literature.

In contrast to the KDEM model, we propose an approach in this paper based on fractional calculations, leading to a process of the dynamics of the evolution of the accumulation and elimination of the contaminant. A process that takes into account the delay in the elimination of the contaminant when the immune system is weakened by a pathogen or a situation of malnutrition in children and pregnant women. However, we propose an FDE that models this delay in elimination. Let us recall for this purpose, that in biology, it was deduced that the membranes of cells of the biological organism have electrical conductance of fractional order [6] and then are classified in the group of models models including noninteger values.

The main contribution of this paper is to propose a dynamic model that takes into account the accumulation and elimination of contaminants from the human body of subjects living with comorbidities. In this context, the exposure process is described by a random variable X, which represents the quantity of the contaminant ingested during a short period.

The rest of the paper is organized as follows: In section 2, we will first lay the mathematical foundations to build the stochastic process that will be used to dynamically model the evolution of the accumulation and elimination of a contaminant when the human organ system fails. Section 3 proposes a stochastic model called the fractional differential equation (FDE) of order 1/2 whose particularity is to take into account the delay in the evolution of a phenomenon. Moreover, we make a brief comparison of the existing classical KDEM model that we implemented with the dioxin dietary intake data and the first model that we elaborated with the failed immune systems.

2. Materials and Methods

In this section, we summarise the fractional integration of Riemann–Liouville and the Adomian decomposition method. These mathematical tools turn out to be very important for our study.

2.1. An Overview of Fractional Integral of Riemann–Liouville

Successive integer-order derivatives have been used to model many life phenomena, including exposure to food risks. Here, our approach is based on the fractional integral of Riemann–Liouville (see [7, 8] for more details).

Let $\alpha \in {ℝ}_{+}^{\ast },a\in ℝ$ and f a locally integrable function defined on $a,+\infty$, the integral of the order $\alpha$ of $f$ of lower bound a is defined for all $t\ge \alpha ,$ such as$$\begin{array}{}\text{(4)}& {}_a{{\displaystyle I}}_t^{\alpha }ft=\frac{1}{\Gamma \alpha }{\displaystyle {\int }_a^t{t-\tau }^{\alpha -1}}f\tau \mathrm{d}\tau ,\end{array}$$while $\mathrm{\Gamma }x={\displaystyle {\int }_0^{+\infty }{t}^{x-1}}{e}^{-t}\mathrm{d}t$ is the well-known Euler gamma function. This integral is convergent, for all $x>0.$ Note that the Beta function is related to the Gamma one and plays an important role in fractional calculations which is defined by [9] such as, for all x, y $>rbin$ 0,$$\begin{array}{}\text{(5)}& Bx,y={\displaystyle {\int }_0^1{t}^{x-1}}{1-t}^{y-1}dt=\frac{\Gamma x\Gamma y}{\Gamma x+y}.\end{array}$$

Furthermore, the fractional left derivative of order $\alpha$ over a given interval $a,b$ of Riemann–Liouville is such as, for all $\alpha >0,x\in a,b$ and n =  $\alpha +1$:$$\begin{array}{}\text{(6)}& {\mathcal{D}}_{a^{+}}^{\alpha }fx=\frac{1}{\Gamma n-\alpha }\frac{d^n}{d{x}^n}{\displaystyle {\int }_a^x{x-t}^{n-\alpha -1}}ftdt,\end{array}$$where $\alpha$ is the integer part of $\alpha$.

2.2. An Overview of the Adomian Decomposition Method

Suppose that we need to solve a functional equation of the form,$$\begin{array}{}\text{(7)}& Au=f;whereA:H⟶H,\end{array}$$is a nonlinear operator defined on a real Hilbert space H, u is the unknown function defined in H, and f is a function given in H. First, we decompose the operator A into$$\begin{array}{}\text{(8)}& A=L+R+N,\end{array}$$where L is an operator inversible in the sense of Adomian, R is the rest of the linear part, and N is the nonlinear part. So, (7) provides the relation as follows:$$\begin{array}{}\text{(9)}& Lu+Ru+Nu=f.\end{array}$$

The operator L is inversible, and we obtain the Adomian canonical form [10].$$\begin{array}{}\text{(10)}& u=\theta +L^{-1}f-L^{-1}Ru-L^{-1}Nu,\end{array}$$where $\theta$ is a constant such that $L\theta =0.$ In particular, the solution of (7) is sought in the form of a series, and we decompose the nonlinear part into$$\begin{array}{}\text{(11)}& Nu={\displaystyle \sum _{n=0}^{+\infty }{A}_k{u}_0,\dots ,u_k}.\end{array}$$

Furthermore, one has$$\begin{array}{}\text{(12)}& {\displaystyle \sum _{n=0}^{+\infty }{u}_n}=\theta +L^{-1}f-L^{-1}R{\displaystyle \sum _{n=0}^{+\infty }{u}_n}-L^{-1}N{\displaystyle \sum _{n=0}^{+\infty }{u}_n}.\end{array}$$

So, we obtain the terms the following numerical series, provideding convergence${\displaystyle {\sum }_{n=0}^{+\infty }{u}_n}$; for all $n\ge 0$, see [11].$$\begin{array}{}\text{(13)}& \begin{array}{c}{u}_0=\theta +L^{-1}f\hfill \\ u_{n+1}=-L^{-1}R{u}_n-L^{-1}N{u}_n\hfill \end{array}.\end{array}$$

The nonlinear part is obtained from the Adomian polynomials defined as follows:$$\begin{array}{}\text{(14)}& A_n=\frac{1}{n!}{\frac{d^n}{d{\lambda }^n}N{\displaystyle \sum _{i=1}^{+\infty }{\lambda }^i}{u}_i}_{\lambda =0}.\end{array}$$

In practice, the following formula is used to calculate polynomials (see [12] for more details). For all n$\ge 0$, one has$$\begin{array}{}\text{(15)}& \begin{array}{c}{A}_0=N{u}_0,\hfill \\ A_{n+1}=\frac{1}{n+1}{\displaystyle \sum _{k=0}^{n}k+1}{u}_{k+1}\frac{\partial A_n}{\partial u_k}.\hfill \end{array}\end{array}$$

The following section deals with the main of our paper results:

3. Main Results

In this paper, we propose an alternative method to the existing ones by the noninteger order Riemann–Liouville derivation, an integration that is the inverse operation of the derivation, which effectively accounts for the delay in the elimination of the contaminant from the body by the immune system during the evolution of the phenomenon. Before building the model, let us make the following assumptions:

  (C1): The subjects living with comorbidities have a less efficient immune system and therefore fail in contrast to immunocompetent subjects.

  (C2): The amount of food toxins in the body today depends on the amount accumulated yesterday. So, the amount of contaminations depends on the time.

Let us notice that when the human organism is immunocompromised, the rate of elimination of the contaminant from the body by the immune system is slow compared to an immunocompetent organism. Hence, there is a need for a fractional differential equation that takes into account the delay in the evolution of the phenomenon. Let us begin with this example where we present the difference between the noninteger derivation and the fractional derivation in the sense of Riemann–Liouville of function is that the noninteger one undergoes a delay in its evolution.$$\begin{array}{}\text{(16)}& D_{0^{+}}^{1/2}x=\frac{1}{\Gamma 1/2}\frac{d}{dx}{\displaystyle {\int }_0^x{x-t}^{-1/2}}tdt,\end{array}$$where $D_{0^{+}}^{1/2}x$ is the variation in the amount of contaminants in the body. Now, by using the well-known relationship, the change of variable $dt=xds,s\in \mathrm{0,1}$.

It comes that$$\begin{array}{}\text{(17)}& D_{0^{+}}^{1/2}x=\frac{1}{\Gamma 1/2}\frac{d}{dx}{\displaystyle {\int }_0^1{x-sx}^{-1/2}}sx.xds,\end{array}$$ furthermore, after a simple calculation, one obtains$$\begin{array}{}\text{(18)}& D_{0^{+}}^{1/2}x=\frac{3}{2}{x}^{1/2}\frac{1}{\Gamma 1/2}\frac{\Gamma 1/2\Gamma 2}{\Gamma 5/2}=\frac{2\sqrt{x}}{\sqrt{\pi }}.\end{array}$$

Remark 1.

A relatively delayed stochastic process can be modeled by a fractional differential equation.

3.1. Univariate Stochastic Process

Note that when a system is disturbed or delayed, it allows a slowed motion; hence, there is a need to model its trajectory by fractional derivation. So, subjects living with comorbidities develop pathology more quickly following exposure to food toxins. This is due to the slow elimination of toxins from their immune systems. The delay in the elimination of contaminants ingested during dietary intake in comorbid individuals can be wellmodeled by a fractional differential equation of order 1/2. Indeed, the linear one-compartment pharmacokinetic model long used by toxicologists is that of an ordinary differential equation of order 1.

Proposition 1.

Let A be the initial body burden of the contaminant at date $T_0=0$. Between two consecutive intakes, the amount of contaminant in the immunodeficient organism decreases with time according to the following fractional differential equation of noninteger order:$$\begin{array}{}\text{(19)}& \begin{array}{c}{D}^{1/2}Xt=-\theta Xt,\hfill \\ X0=A.\hfill \end{array}\end{array}$$

The exposure computed immediately after the i-th food intake is such as$$\begin{array}{}\text{(20)}& Xt=A{e}^{{\theta }^{2}t}-\frac{4\theta \sqrt{t}}{\sqrt{\pi }}{\displaystyle \sum _{n=0}^{+\infty }\frac{4^{n}n+1!{{\theta }^{2}t}^n}{2n+2!}},\end{array}$$where $Xt$ is the amount of the contaminant in the body at a given time t, $\theta >0$ is the removal rate of the contaminant between the dates $T_i$ et $T_{i+1},$ where $\theta =ln2/DV,$ DV is the biological half-life of the contaminant, and X(0) is the initial body burden of the contaminant at $T_0=0.$

The human organism has a system of accumulation in contaminants, and these last ones are eliminated in a progressive way and also depend on the frequency of the food intake.

Proof.

We propose a proof based on fractional Riemann–Liouville integration and the Adomian decomposition method.

Consider the following equation:$$\begin{array}{}\text{(21)}& \begin{array}{c}{D}^{1/2}Xt+\theta Xt=0\hfill \\ X0=A\hfill \end{array}.\end{array}$$

Then, it follows that$$\begin{array}{}\text{(22)}& \begin{array}{c}LX=D_{0^{+}}^{1/2}X\\ RX=\theta Xt\\ g=0\end{array}\text{where\hspace{0.17em}}{L}^{-1}=I_{0^{+}}^{1/2}.\end{array}$$

By using (21), it comes that$$\begin{array}{}\text{(23)}& LX+RX=g.\end{array}$$

By using $L^{-1}$ to (23), one has$$\begin{array}{}\text{(24)}& L^{-1}LX=I_{0^{+}}^{1/2}{D}_{0^{+}}^{1/2}Xt=Xt-X0.\end{array}$$

It comes that$$\begin{array}{}\text{(25)}& Xt-X0+L^{-1}RX=L^{-1}g.\end{array}$$

That allows us to consider the form:$$\begin{array}{}\text{(26)}& X={\displaystyle \sum _{n=0}^{+\infty }{X}_n}.\end{array}$$

From (25) and (26), it yields that$$\begin{array}{}\text{(27)}& {\displaystyle \sum _{n=0}^{+\infty }{X}_n}t=X0+L^{-1}g-{\displaystyle \sum _{n=0}^{+\infty }{L}^{-1}R{X}_{n}t},\end{array}$$giving the following system:$$\begin{array}{}\text{(28)}& \begin{array}{c}{X}_{0}t=X0+L^{-1}g,\\ X_{n+1}t=-L^{-1}R{X}_{n}tn\ge 0.\end{array}\end{array}$$

Calculation of the terms of $X_{n}t.$

At order 0, we have $X_{0}t$ such as$$\begin{array}{}\text{(29)}& X_{0}t=X0+L^{-1}g=A.\end{array}$$

In order 1, $X_{1}t$, using (29) one has$$\begin{array}{}\text{(30)}& X_{1}t==-I_{0^{+}}^{1/2}\theta X_{0}t=-I_{0^{+}}^{1/2}\theta A=-\frac{2\theta A{t}^{1/2}}{\sqrt{\pi }}.\end{array}$$By using (30), one has$$\begin{array}{}\text{(31)}& X_{2}t=-L^{-1}R{X}_{1}t=-I_{0^{+}}^{1/2}\theta \frac{-2\theta A{t}^{1/2}}{\sqrt{\pi }}=\frac{2{\theta }^{2}A{I}_{0^{+}}^{1/2}{t}^{1/2}}{\sqrt{\pi }}.\end{array}$$

From (31) and (5), it comes that$$\begin{array}{}\text{(32)}& I_{0^{+}}^{1/2}{t}^{1/2}=\frac{1}{\Gamma 1/2}{\displaystyle {\int }_0^1{t-st}^{-1/2}}{st}^{1/2}tds=\frac{t}{\Gamma 1/2}{\displaystyle {\int }_0^1{1-s}^{1/2-1}}{s}^{1/2}ds=\frac{t\sqrt{\pi }}{2}.\end{array}$$

By introducing (32) into (31), it comes that$$\begin{array}{}\text{(33)}& X_{2}t=\frac{2{\theta }^{2}A}{\sqrt{\pi }}\frac{t\sqrt{\pi }}{2}={\theta }^{2}At.\end{array}$$

For $X_{3}t$, one has$$\begin{array}{}\text{(34)}& X_{3}t=-L^{-1}R{X}_{2}t.\end{array}$$

By replacing (33) in the above formula, it comes that$$\begin{array}{}\text{(35)}& X_{3}t=-I_{0^{+}}^{1/2}\theta {\theta }^{2}At=-{\theta }^{3}A{I}_{0^{+}}^{1/2}t=\frac{t^{1/2}}{\Gamma 1/2}{\displaystyle {\int }_0^1{1-s}^{-1/2}}sds=\frac{-4{\theta }^{3}A{t}^{3/2}}{3\sqrt{\pi }}.\end{array}$$

In an analogous way, one obtains successively the terms$$\begin{array}{}\text{(36)}& \begin{array}{c}{X}_{0}t=A\hfill \\ X_{1}t=-\frac{2\theta A{t}^{1/2}}{\sqrt{\pi }}\hfill \\ X_{2}t={\theta }^{2}At\hfill \\ X_{3}t=\frac{-4{\theta }^{3}A{t}^{3/2}}{3\sqrt{\pi }}\hfill \\ X_{4}t=\frac{{\theta }^{4}A{t}^2}{2}\hfill \\ X_{5}t=-\frac{8{\theta }^{5}A{t}^{5/2}}{15\sqrt{\pi }}\hfill \\ X_{6}t=\frac{{\theta }^{6}A{t}^3}{6}\hfill \\ X_{7}t=-\frac{16{\theta }^{7}A{t}^{7/2}}{105\sqrt{\pi }}\hfill \\ ⋮\hfill \end{array}.\end{array}$$

However, in practice, all terms of the series cannot be determined, so we use an approximation of the solution from the truncated series as follows:$$\begin{array}{}\text{(37)}& Xt={\displaystyle \sum _{n=0}^{+\infty }{X}_{n}t},\end{array}$$which provides the following equivalences:$$\begin{array}{c}\text{(38)}& Xt\simeq X_{0}t+X_{1}t+X_{2}t+X_{3}t+X_{4}t+X_{5}t+X_{6}t+X_{7}t\cdots \simeq A-\frac{2\theta A{t}^{1/2}}{\sqrt{\pi }}+{\theta }^{2}At-\frac{4{\theta }^{3}A{t}^{3/2}}{3\sqrt{\pi }}+\frac{{\theta }^{4}A{t}^2}{2}-\frac{8{\theta }^{5}A{t}^{5/2}}{15\sqrt{\pi }}+\frac{{\theta }^{6}A{t}^3}{6}-\frac{16{\theta }^{7}A{t}^{7/2}}{105\sqrt{\pi }}+\cdots \\ \simeq A1+{\theta }^{2}t+\frac{{\theta }^4{t}^2}{2!}+\frac{{\theta }^6{t}^3}{3!}+\cdots +\\ +\frac{A}{\sqrt{\pi }}-\frac{\theta t^{1/2}}{1/2}-\frac{{\theta }^3{t}^{3/2}}{3/2\times 1/2}-\frac{{\theta }^5{t}^{5/2}}{5/2\times 3/2\times 1/2}-\frac{{\theta }^7{t}^{7/2}}{7/2\times 5/2\times 3/2\times 1/2}+\cdots .\end{array}$$

In the same way, it follows that$$\begin{array}{c}\text{(39)}& Xt\simeq A{\displaystyle \sum _{n=0}^{+\infty }\frac{{{\theta }^{2}t}^n}{n!}}-A{\displaystyle \sum _{n=0}^{+\infty }\frac{{{\theta }^{2}t}^{1/2+n}}{\Gamma 1/2+n+1}}\simeq A{e}^{{\theta }^{2}t}-A{\displaystyle \sum _{n=0}^{+\infty }\frac{{{\theta }^{2}t}^{1/2+n}}{\Gamma 1/2+n+1}}\\ \simeq A{e}^{{\theta }^{2}t}-A\frac{4\sqrt{{\theta }^{2}t}}{\sqrt{\pi }}{\displaystyle \sum _{n=0}^{+\infty }\frac{4^{n}n+1!{{\theta }^{2}t}^n}{2n+2!}}\\ \simeq A{e}^{{\theta }^{2}t}-A\frac{4\theta \sqrt{t}}{\sqrt{\pi }}{\displaystyle \sum _{n=0}^{+\infty }\frac{4^{n}n+1!{{\theta }^{2}t}^n}{2n+2!}}\cdot \end{array}$$

Finally, one obtains$$\begin{array}{}\text{(40)}& Xt\simeq A{e}^{{\theta }^{2}t}-\frac{4\theta \sqrt{t}}{\sqrt{\pi }}{\displaystyle \sum _{n=0}^{+\infty }\frac{4^{n}n+1!{{\theta }^{2}t}^n}{2n+2!}}.\end{array}$$

Note that $Xt$ represents the behavior of the quantity of contaminated intake at each consumption in the human body through time for subjects living with comorbidities.

Proposition 2.

Under (C1), the trajectory of the accumulated exposure is defined as follows:$$\begin{array}{}\text{(41)}& X_{i+1}t=X_i{e}^{{\theta }^{2}t}-\frac{4\theta \sqrt{t}}{\sqrt{\pi }}{\displaystyle \sum _{n=0}^{+\infty }\frac{4^{n}n+1!{{\theta }^{2}t}^n}{2n+2!}}+E_{i+1}.\end{array}$$$E_{i+1}$ is the calculated exposure at time $t_{i+1}.$

Proof.

According to proposition 1, the first term of (41) is proved. We obtain the second term proposed in the work of Bertail in [4, 5], i.e., the exposure over a given period is none other than the sum of the consumption of the contaminated products, weighted by the contamination rates associated with each of the products of doubtful quality. By noting P as the number of foodstuffs carrying the contamination $C_P$, the consumption of any individual of body weight $w$, and $Q_P$ is the contamination rate in $\mathrm{\mu g}/\text{kg}$ or pg/kg of each of these foodstuffs. This leads to the result of the random exposure of an individual to the contaminant of interest as follows (3) concluding the proof.

3.2. Conditional Dependence of the Multivariate Stochastic Process

Let us make some additional assumptions.

  (C3): Let us define the health risk arising from the accumulation and interaction between the different contaminants in the body. For instance, this involves assessing the risk due to the accumulation of three substances in the body. It should be noted that studies have shown that there is a dependency between the trajectories of contaminants in the body.

  (C4): The body burden of each contaminant during the month depends on the body burden of that contaminant in the previous month and the new exposure.

Based on the fact that contaminants are stored and eliminated by the same organ, we make the following dependency hypothesis. The above assumptions lead to the following proposition:

Proposition 3.

Let $X_i=X_{i1},X_{i2},X_{i3},i=1,\dots ,n$ be the vector body burden of each contaminant during $i^{th}$ week or month. Then, the evolution of the quantity of the contaminant over time is obtained by the following relation:$$\begin{array}{}\text{(42)}& X_i=X_{i-1}\psi +E_i,\end{array}$$with $\psi ={\mathrm{ℳ}}_{3ℝ},$ the matrix of transition probabilities is given by$$\begin{array}{}\text{(43)}& \begin{array}{ccc}{e}^{{\theta }_1^{2}t}-\frac{4{\theta }_1\sqrt{t}}{\sqrt{\pi }}{\displaystyle \sum _{n=0}^{+\infty }\frac{4^{n}n+1!{{\theta }_1^{2}t}^n}{2n+2!}}& a& b\\ c& e^{{\theta }_2^{2}t}-\frac{4{\theta }_2\sqrt{t}}{\sqrt{\pi }}{\displaystyle \sum _{n=0}^{+\infty }\frac{4^{n}n+1!{{\theta }_2^{2}t}^n}{2n+2!}}& d\\ e& f& e^{{\theta }_3^{2}t}-\frac{4{\theta }_3\sqrt{t}}{\sqrt{\pi }}{\displaystyle \sum _{n=0}^{+\infty }\frac{4^{n}n+1!{{\theta }_3^{2}t}^n}{2n+2!}}\end{array},\end{array}$$where ${\theta }_1,{\theta }_2,{\theta }_3$ denote the removal rates of the contaminant in the body and $E_i=E_{i1},E_{i2},E_{i3}$ the vectors of the trajectories of the exposures to the contaminants during the $i^{th}$ week or month. The matrix $\psi$ contains the contaminant removal coefficients.

Proof.

The proof of the previous proposition based on assumptions (C3) and (C4) yield the desired result.

It can be seen easily that$$\begin{array}{}\text{(44)}& \Rightarrow \begin{array}{c}{X}_{i1}=X_{i-\mathrm{1,1}}\ast e^{{\theta }_1^{2}t}-\frac{4{\theta }_1\sqrt{t}}{\sqrt{\pi }}{\displaystyle \sum _{n=0}^{+\infty }\frac{4^{n}n+1!{{\theta }_1^{2}t}^n}{2n+2!}}+c\ast X_{i-\mathrm{1,2}}+e\ast X_{i-\mathrm{1,3}}+E_{i1}\\ X_{i2}=a\ast X_{i-\mathrm{1,1}}+X_{i-\mathrm{1,2}}\ast e^{{\theta }_2^{2}t}-\frac{4{\theta }_2\sqrt{t}}{\sqrt{\pi }}{\displaystyle \sum _{n=0}^{+\infty }\frac{4^{n}n+1!{{\theta }_2^{2}t}^n}{2n+2!}}+f\ast X_{i-\mathrm{1,3}}+E_{i2}\\ X_{i3}=b\ast X_{i-\mathrm{1,1}}+d\ast X_{i-\mathrm{1,2}}+X_{i-\mathrm{1,3}}\ast e^{{\theta }_3^{2}t}-\frac{4{\theta }_3\sqrt{t}}{\sqrt{\pi }}{\displaystyle \sum _{n=0}^{+\infty }\frac{4^{n}n+1!{{\theta }_3^{2}t}^n}{2n+2!}}+E_{i3}\end{array}.\end{array}$$

Neglecting the impact of the removal of one contaminant on the other, we have trajectory (41) for each contaminant considered, i.e., the coefficients a, b, c, and d are zero. However, when the contaminants are stored in the same organ, then the elimination becomes slower because of the dependency of the elimination. Assumptions will be made to determine the reals a, b, c, d, e, and f of the matrix by an Archimedean copula dependence investigation in [13, 14] and estimate the probability that the amount of contaminant 1 in the body exceeds a threshold u1 knowing that contaminant 2 has exceeded the threshold u2. Unfortunately, this will not be examined in this work.

3.3. Generalized Fractional Differential Equation

Assume once again that the manifestation of a pathology differs from one person to another. Based on this assumption, we propose a generalized fractional differential equation of order $\alpha$ which takes into account the delay in the evolution of a phenomenon.

Proposition 4.

The following system provides a fractional differential equation:$$\begin{array}{}\text{(45)}& \begin{array}{c}{D}^{\alpha }Xt+\theta Xt=0\hfill \\ X0=\kappa \hfill \end{array},\end{array}$$such that t$>$ 0; $0<\alpha \le 1$; modeling the evolution through time of the dynamics of accumulation and elimination of the contaminant in the human organism with a deficient immune system during consecutive intakes of the contaminated food. Then, the solution is given by$$\begin{array}{}\text{(46)}& Xt={\displaystyle \sum _{n=0}^{+\infty }\frac{{-1}^n{\theta }^n\kappa t^{n\alpha }}{\Gamma n\alpha +1}},\end{array}$$It is a stochastic time dependent process that quantifies the exposure to toxins of subjects living with comorbidities.

Proof.

The proof requires frequent uses of parts section 2.1 and section 2.2 of section 2. Indeed, it is easy to show that$$\begin{array}{}\text{(47)}& X_{0}t=X0=\kappa .\end{array}$$

In order 1, the quantity $X_{1}t,$ gives$$\begin{array}{}\text{(48)}& X_{1}t=-L^{-1}R{X}_{0}t=-\frac{\theta \kappa t^{\alpha }}{\Gamma \alpha +1}.\end{array}$$

By introducing (48) into the following relation, one has$$\begin{array}{}\text{(49)}& X_{2}t=-L^{-1}R{X}_{1}t=\frac{{\theta }^2\kappa I_{0^{+}}^{\alpha }{t}^{\alpha }}{\Gamma \alpha +1}.\end{array}$$

So, we get the following results:$$\begin{array}{}\text{(50)}& I_{0^{+}}^{\alpha }{t}^{\alpha }=\frac{1}{\Gamma \alpha }{\displaystyle {\int }_0^1{t-st}^{\alpha -1}}{st}^{\alpha }tds=\frac{t^{2\alpha }}{\Gamma \alpha }\beta \alpha ,\alpha +1=\frac{t^{2\alpha }\Gamma \alpha +1}{\Gamma 2\alpha +1}.\end{array}$$

Furthermore, substituting (50) into (49) one has$$\begin{array}{}\text{(51)}& X_{2}t=\frac{{\theta }^2\kappa }{\Gamma \alpha +1}\times \frac{t^{2\alpha }\Gamma \alpha +1}{\Gamma 2\alpha +1}=\frac{{\theta }^2\kappa t^{2\alpha }}{\Gamma 2\alpha +1}.\end{array}$$

For $X_{3}t,$ we have$$\begin{array}{}\text{(52)}& X_{3}t=-L^{-1}R{X}_{2}t=-I_{0^{+}}^{\alpha }\frac{{\theta }^3\kappa t^{2\alpha }}{\Gamma 2\alpha +1}=\frac{-{\theta }^3\kappa }{\Gamma 2\alpha +1}{I}_{0^{+}}^{\alpha }{t}^{2\alpha }.\end{array}$$

Applying the same operator, it comes that$$\begin{array}{}\text{(53)}& I_{0^{+}}^{\alpha }{t}^{2\alpha }=\frac{1}{\Gamma \alpha }{\displaystyle {\int }_0^1{t}^{3\alpha }}{1-s}^{\alpha -1}{s}^{2\alpha }ds=\frac{t^{3\alpha }}{\Gamma \alpha }\beta \alpha ,2\alpha +1.\end{array}$$

And finally, it comes that$$\begin{array}{}\text{(54)}& X_{3}t=\frac{-{\theta }^3\kappa }{\Gamma 2\alpha +1}\frac{t^{3\alpha }\Gamma 2\alpha +1}{\Gamma 3\alpha +1}=\frac{-{\theta }^3\kappa t^{3\alpha }}{\Gamma 3\alpha +1}.\end{array}$$

By conjecture, we have the solution to the problem (21) as follows:$$\begin{array}{}\text{(55)}& X_{n}t=\frac{{-1}^n{\theta }^n\kappa t^{n\alpha }}{\Gamma n\alpha +1}.\end{array}$$

According to (10), we finally obtain the following result:$$\begin{array}{}\text{(56)}& Xt={\displaystyle \sum _{n=0}^{+\infty }{X}_n}t={\displaystyle \sum _{n=0}^{+\infty }\frac{{-1}^n{\theta }^n\kappa t^{n\alpha }}{\Gamma n\alpha +1}}.\end{array}$$

Our result is proved if we show by recurrence that$$\begin{array}{}\text{(57)}& Xt=\frac{{-1}^n{\theta }^n\kappa t^{n\alpha }}{\Gamma n\alpha +1}.\end{array}$$

For the order 0, $n=0,Xt=X_{0}t=\kappa$.

The property is true at order $n=0$; now let us suppose it is true at any order $n$ and show that the property is true at order $n+1$, i.e.,$$\begin{array}{}\text{(58)}& X_{n+1}t=\frac{{-1}^{n+1}{\theta }^{n+1}\kappa t^{n+1\alpha }}{\Gamma n+1\alpha +1}.\end{array}$$

Using the same approach, we obtain the results successively as follows:$$\begin{array}{}\text{(59)}& X_{n+1}t=-L^{-1}R{X}_{n}t=-I_{0^{+}}^{\alpha }\theta \frac{{-1}^n{\theta }^n\kappa t^{n\alpha }}{\Gamma n\alpha +1},\end{array}$$which gives,$$\begin{array}{}\text{(60)}& X_{n+1}t==\frac{{-1}^{n+1}{\theta }^{n+1}\kappa I_{0^{+}}^{\alpha }{t}^{n\alpha }}{\Gamma n\alpha +1}.\end{array}$$

Thus,$$\begin{array}{}\text{(61)}& I_{0^{+}}^{\alpha }{t}^{n\alpha }=\frac{1}{\Gamma \alpha }{\displaystyle {\int }_0^1{t-st}^{\alpha -1}}{st}^{n\alpha }tds=\frac{t^{n+1}}{\Gamma \alpha }{\displaystyle {\int }_0^1{1-s}^{\alpha -1}}{s}^{n\alpha }ds.\end{array}$$

This is equivalent to$$\begin{array}{}\text{(62)}& I_{0^{+}}^{\alpha }{t}^{n\alpha }=\frac{1}{\Gamma \alpha }\beta \alpha ,n\alpha +1=\frac{t^{n+1}\alpha \Gamma n\alpha +1}{\Gamma n+1\alpha +1},\end{array}$$and consequently,$$\begin{array}{}\text{(63)}& X_{n+1}t=\frac{{-1}^{n+1}{\theta }^{n+1}\kappa }{\Gamma n\alpha +1}\frac{t^{n+1\alpha }\Gamma n\alpha +1}{\Gamma n+1\alpha +1}=\frac{{-1}^{n+1}{\theta }^{n+1}\kappa t^{n+1\alpha }}{\Gamma n+1\alpha +1}.\end{array}$$

Thus, the property is also true to order $n+1$. So, it is true to any order $n$, i.e.,$$\begin{array}{}\text{(64)}& X_{n}t=\frac{{-1}^n{\theta }^n\kappa t^{n\alpha }}{\Gamma n\alpha +1}.\end{array}$$

Finally, one has$$\begin{array}{}\text{(65)}& Xt={\displaystyle \sum _{n=0}^{+\infty }{X}_n}t={\displaystyle \sum _{n=0}^{+\infty }\frac{{-1}^n{\theta }^n\kappa t^{n\alpha }}{\Gamma n\alpha +1}}.\end{array}$$

Remark 2.

By considering $\alpha =1/2$, we obtain the formula by (28).

The paper is closed with a comparison and simulation regarding the new model FDE and the existing one KDEM extensively commented on in the introduction.

3.4. Comparative Study and Simulation

The main part of this section is to extend previous theoretical investigation results (15) graphically. The comparison of the two models, the existing one and the one we propose, is obtained by application to dioxins. The elimination rate is $\theta$  = 0.006418 and assume that the initial body burden at $T_0=0$ is $X0=15\hspace{0.17em}\text{pg}/\text{kg}$.

Table 1 shows that our proposed process results in a slightly higher amount of the contaminant over time compared to the existing model. This proves that this model is well adapted to failing immune systems. Thus, it is well suited to model the dynamics of contaminant evolution in the body of people living with comorbidities, fragile children, and pregnant women. Since the KDEM model is adapted to the nonfailing immune system, we therefore implemented it for a graphical view of the evolution of the contaminant dynamics in the immunocompetent organism.

Table 1

Evolution of the quantity of the contaminant over time.

Time (months)	$X_{i+1}t=X_i.e^{-\theta t}$	$X_{i+1}t=X_i{e}^{{\theta }^{2}t}-4\theta \sqrt{t}/\sqrt{\pi }{\displaystyle {\sum }_{n=0}^{+\infty }{4}^{n}n+1!{{\theta }^{2}t}^n/2n+2!}$
0	15	15
1	14,99903733	14,99891376
2	14,99711218	14,99737774
3	14,99422491	14,99549674

The graph above Figure 1 shows the evolution of the dynamics of dioxins once introduced into the human organism. The value of the elimination rate (0.006418) is obtained by the formula $\theta =\text{ln}2/DV$ (DV the biological half-life of dioxin is 9 years or 108 months). Given that the tolerable threshold of dioxin by the human organism is 70pg/kg of body weight, we assume however the initial body load at 10 pg/kg of body weight. We note that in this process the time taken to exceed the tolerable threshold is long, in contrast to the trajectory of contaminants in the body of people living with a comorbidity, where this threshold is exceeded more quickly after contamination Figure 2. We recall that this exceeding can prove to be dangerous for the organism.

[figure(s) omitted; refer to PDF]

This graphical test shows that organisms living with comorbidities, people with a weak immune system, namely the child not breastfed or malnourished, pregnant women, all these people, the system eliminates the contaminant slowly, which exposes the danger of developing a pathology more quickly than immunocompetent subjects.

4. Conclusion and Discussions

We have developed a model that takes into account a weak immune system for modeling the risks associated with poor nutrition. This allows medical services to properly assess and prevent the danger of exposure to food toxins and improve the life expectancy of these people already weakened by pathologies such as asthma, diabetes, cardiovascular diseases, AIDS, especially malnourished children, people with disabilities, people who are not well nourished, and people of advanced age and pregnant women. We have shown that the toxin values in the organism successively obtained by the FDE model are slightly higher than the amount obtained by the existing model. Naturally, it is clear that our model is better than the KDEM model, but this confirms the hypothesis (C1).

This fractional differential equation model that we propose in this paper to quantify dietary risk exposure marks’ progress in the analysis and search for solutions for the prevention of the above-mentioned diseases. In our next research, we will take into account in this new FDE model, the multiple contaminations to which the organism is exposed daily, especially the application to contaminants found in food consumption in sub-Saharan African countries.

Acknowledgments

The authors would like to thank the Academic Editor of Hindawi Journal and all the reviewers of this article for their comments and suggestions. Thank you for all that you do to advance science. Université Thomas Sankara, Burkina Faso.

References

[1] Ph Verger, J. Tressou, S. Clémençon, "Integration of time as a description parameter in risk characterisation: application to methyl mercury," Regulatory Toxicology and Pharmacology, vol. 49, pp. 25-30, DOI: 10.1016/j.yrtph.2007.04.010, 2007.

[2] C. Tillier, Théorie de la ruine et risque alimentaire, 2014. http://papersjds14.sfds.asso.fr/submission_254.pdf

[3] J. Tressou, A. Crépet, P. Bertail, M. H. Feinberg, J. Ch Leblanc, "Probabilistic exposure assessment to food chemicals based on extreme value theory. Application to heavy metals from fish and sea products," Food and Chemical Toxicology, vol. 42, pp. 1349-1358, DOI: 10.1016/j.fct.2004.03.016, 2004.

[4] S. Clémençon, J. Tressou, "Exposition aux risques alimentaires et processus stochastiques: le cas des contaminants chimiques," Journal de la Société Française de Statistique, vol. 15, 2009. http://www.numdam.org/item/JSFS_2009__150_1_3_0/

[5] F. Debes, E. Budtz-Jorgensen, P. Weihe, R. F. White, P. Grandjean, "Impact of prenatal methylmercury exposure on neurobehavioral function at age 14 years," Neurotoxicology and Teratology, vol. 28, pp. 536-547, DOI: 10.1016/j.ntt.2006.02.005, 2006.

[6] G. Chen, E. Friedman, "An RLC interconnect model based on Fourier analysis," IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 24 no. 2, pp. 170-183, DOI: 10.1109/TCAD.2004.841065, 2005.

[7] J. M. Loufouilou, G. Bissanga, B. Francis, Y. Pare, "Application of Adomian decomposition method to solving some kinds of partial differential equations and system of partial differential equations," International Journal of Advanced Mathematical Sciences, Citeseer, vol. 4, pp. 190-198, DOI: 10.1.1.1012.3799, 2013.

[8] T. Abdeljawad, M. Grossman, "On geometric fractional calculus," J. Semigroup Theory Appl., vol. 2016, 2016. http://www.scik.org/index.php/jsta/article/view/2546

[9] K. Fridjet, Dérivée et intégrale fractionnaire au sens de Riemann-Liouville, 2018. http://dspace.univ-eloued.dz/handle/123456789/517

[10] J. B. Yindoula, Y. Pare, G. Bissanga, F. Bassono, B. Some, "Application of the Adomian decomposition method and Laplace transform method to solving the convection diffusion-dissipation equation," International Journal of Applied Mathematical Research, vol. 3 no. 1,DOI: 10.1.1.970.3228, 2013.

[11] B. Béyi, J. Bonazebi-Yindoula, L. Somé, G. Bissanga, "Application of the laplace-adomian method and the SBA method to solving the partial differential and integro-differential equations," Journal of Mathematical Sciences: Advances and Applications,DOI: 10.18642/jmsaa_7100122118, .

[12] S. A. El-Wakil, M. A. Abdou, "New applications of Adomian decomposition method, Chaos," Solitons & Fractals, vol. 33, pp. 513-522, DOI: 10.1016/j.chaos.2005.12.037, 2007.

[13] D. Barro, M. Diallo, R. G. Bagré, "Spatial tail dependence and survival stability in a class of archimedean copulas," International Journal of Mathematics and Mathematical Sciences, vol. 2016,DOI: 10.1155/2016/8927248, 2016.

[14] S. Dossou-Gbete, B. Some, D. Barro, "Modelling the dependence of parametric bivariate extreme value copulas," Asian Journal of Mathematics & Statistics, vol. 2 no. 3,DOI: 10.3923/ajms.2009.41.54, 2009.