On the Construction of Some Fractional Stochastic Gompertz Models

: The aim of this paper is the construction of stochastic versions for some fractional Gompertz curves. To do this, we ﬁrst study a class of linear fractional-integral stochastic equations, proving existence and uniqueness of a Gaussian solution. Such kinds of equations are then used to construct fractional stochastic Gompertz models. Finally, a new fractional Gompertz model, based on the previous two, is introduced and a stochastic version of it is provided.


Introduction
Fractional calculus is presently applied to a lot of scientific fields. Despite the problem of defining fractional derivatives being quite old (see, for instance, [1,2]), it has mainly been developed in recent times (see [3]). Due to its versatility in describing slower or also different time scales, fractional derivatives and fractional-order differential equations are very often used in applications, so that also different books have been written on the argument (see, for instance, [4][5][6]). The main generalization of the classical Cauchy problems to the fractional order is achieved via the so-called Caputo-fractional derivative, introduced by Michele Caputo in [7]. In such paper, the fractional derivative is used to study the Q-factor of some non-ferromagnetic solids, thus being introduced in an applicative context. From such moment, fractional calculus has been used to address a lot of different models: from epidemics [8] to osmosis [9], from neurophysiology [10] to viscoelasticity [11] and many others [12].
Here we focus on fractional-order population growth models. A first model of population growth can be achieved by modifying the classical Malthus model by introducing a fractional-order derivative in place of the classical one (see [12,13]). As a second step, one could ask for a fractional-order generalization of a Gompertz model. Gompertz model are quite popular growth model. Such models take into account a time-varying birth rate, which describes the fact that a person's resistance to death decreases with age. Such models have been used in particular to model cancer growth, starting from [14] and then used to describe a single species growth (see for instance [15]). For this and other reasons, Gompertz curves have been widely studied. For instance, knowing that some species of cancer evolved following a Gompertz law, optimal control of it has become necessary (see for instance [16]). At the same time, stochastic models became necessary to describe eventual environmental (and thus unpredictable) effects (see [17][18][19][20] and many others).
Concerning fractional-order Gompertz models, the first one has been introduced in [21], but it is not achieved by simply substituting the fractional derivative in place of the classical one. To understand how the fractional-order model is introduced, let us recall that the classical Gompertz curve x(t) can be defined as the solution of the non-linear Cauchy problem where α, β > 0 are dynamical parameters and x 0 > 0 is the initial population density. It is also well known that the solution is given by and then the model admits a carrying capacity If we define y(t) = log x(t) x 0 , the Cauchy problem (1) can be rewritten as In [21], the fractional-order Gompertz model is achieved by substituting the Caputo-fractional derivative in place of the classical one only in the linear equation. In particular, the function y ν (t) is defined as the solution of the fractional Cauchy problem d ν y ν dt ν (t) = α − βy ν (t), y ν (0) = 0, where d ν dt ν is the fractional Caputo derivative of order ν ∈ (0, 1), and then defining the fractional Gompertz curve as x ν (t) = x 0 e y ν (t) . In such case, we have where E ν (t) is the Mittag-Leffler function (defined in Equation (8)).
In [22] another type of fractional-order Gompertz model has been introduced. To understand how such model is defined, let us recall that for the classical model we have y(t) = α β (1 − e −βt ) and then we can rewrite the first equation of Equations (2) as In [22], they use the Caputo-fractional derivative with respect to another function, as defined in [23], to define the improved fractional Gompertz curve x ν (t) the solution of given by In this paper, we aim to define a class of stochastic Gompertz models that generalize the two proposed fractional Gompertz curves. To do this, we first need to investigate some results related to a class of stochastic linear fractional-integral equations, concerning in particular the existence of Gaussian solutions. Such equations generalize the Caputo-fractional stochastic differential equations studied for instance in [24,25].
In particular this approach leads to a method of construction for general fractional growth models with noise that preserves normal or log-normal one-dimensional distributions. The preservation of such laws permits recognition in some macroscopical observable functions (the mean in the normal case and the median in the log-normal case) of the original growth models. Thus, these stochastic models work as a noisy perturbation of the original deterministic ones. This procedure could not be achieved by using the classical tools of fractionalization via time-change (see for instance [26][27][28][29]) for different reasons. For instance, if we apply a time-change to the stochastic Gompertz model, since the stochastic differential equation that drives the model is non-linear, its mean does not solve Equation (1) with the Caputo derivative in place of the classical one. However, such time-changed process can be still seen as the exponential of a time-changed Ornstein-Uhlenbeck process (in the sense of [28]), but, being the latter not a Gaussian one, the time-changed Gompertz model is not log-normal and its median does not coincide with the function x ν given in Equation (4), despite the mean of the time-changed Ornstein-Uhlenbeck process is still solution of Equation (3). Our new approach overtakes such problems, giving then some log-normal or normal processes whose dynamics are given by perturbation of the deterministic ones.
The paper is structured as follows: • In Section 2 we give some basic definitions and preliminaries on fractional calculus; • In Section 3 we study a class of linear fractional-integral stochastic equations: we will need them to define the stochastic models for fractional Gompertz curves. In particular, we focus on existence and almost surely uniqueness of Gaussian solutions. Moreover, since they are obtained via a Picard approximation method, we also give an estimate of the speed of convergence of the method in terms of the distribution of the maximum of the chosen noise.

•
In Section 4 we give some examples on possible choices of noise. In particular in Section 4.4 we show that such fractional-integral equations are indeed a generalization of the fractional stochastic differential equations discussed in [24,25].

•
In Section 5 we use the results from the previous sections to introduce stochastic models for fractional Gompertz curves. In particular in Section 5.1 we give some generalities on the classical stochastic Gompertz model, while in Sections 5.2 and 5.3 we give a stochastic version of the fractional Gompertz curve introduced in [21] and of the improved fractional Gompertz curve introduced in [22]. Finally, in Section 5.4 we construct a new fractional Gompertz model obtained by merging the approach of the previous two models and we describe a stochastic counterpart for it.

Some Preliminaries on Fractional Calculus
Concerning the main properties of fractional integrals and derivatives, we refer to [30]. Let us give the following definition of the fractional-integral. Definition 1. Given ν > 0 the fractional-integral I ν t of order ν is defined as for any suitable function f : [0, +∞) → R.
It is easy to see that for instance, for any f ∈ L ∞ the fractional-integral I ν t f is defined. Moreover, for any ν 1 , ν 2 it holds I ν 1 t I ν 2 t = I ν 1 +ν 2 t . It is also interesting to notice that the fractional-integral I ν t is a convolution operator. Indeed if we define the kernel I ν (t) = t ν−1 Γ(ν) χ [0,+∞) , then, for any function where * is the convolution product and f is extended to the whole real line by setting f (t) = 0 for any t < 0. If ν ∈ (0, 1), then the convolution kernel I ν is singular, but still in L 1 (0, T) for any T > 0. Therefore, while for any ν ∈ [1, +∞), one only needs f to be in L 1 loc , it is not enough if ν ∈ (0, 1). Now we can define the Riemann-Liouville derivative.

Definition 2.
Given ν ∈ (0, 1), the Riemann-Liouville fractional derivative D ν t of order ν is defined as From now on we will denote f (0+) = f (0). This relation lets us also define the Caputo-fractional derivative for any Riemann-Liouville derivable function, hence for a much wider class of functions. Concerning Caputo derivatives, we can define fractional Cauchy problems by using them. Indeed, under suitable assumptions, the fractional Cauchy problem d ν y dt ν (t) = F(t, y(t)), y(0) = y 0 is well-posed. In particular, the relaxation problem d ν y dt ν (t) = ay(t), y(0) = y 0 admits as unique solution the function which is a generalization of the exponential function (observe that if ν = 1, E 1 (t) = e t ).
We need also to introduce fractional calculus with respect to other functions. Riemann-Liouville type fractional derivative of a function with respect of another function were introduced to deal with Leibniz rule and chain rule for fractional derivatives (see, for instance, [31,32]). For this part, we mainly refer to [23]. Let us first give the definition of fractional-integral with respect to another function. Definition 4. Given ν > 0 and an increasing function Ψ ∈ C 1 (0, +∞) such that Ψ (t) = 0 for any t > 0, the fractional-integral I ν,Ψ t with respect to Ψ is given by for any suitable function f : [0, +∞) → R.
Observe that if Ψ(t) = t, we achieve the classical fractional-integral. Let us now define the Riemann-Liouville type fractional derivative.
Observe that for Ψ(t) = t, we achieve the classical Riemann-Liouville fractional derivative. Moreover, we have in this case Let us also give the definition of the Caputo type fractional derivative.
In [23] (Theorem 3) the following relation is shown Using this relation, one can extend the definition of Caputo-fractional derivative of a function with respect to another function to the whole class of the Riemann-Liouville derivable (with respect to Ψ) functions. Moreover, under suitable assumptions, the following fractional Cauchy problem In the spirit of [31,32], let us show a chain rule for Caputo-fractional derivatives of a function with respect to another function. Proposition 1. Let g be a Caputo-derivable function and Ψ be an increasing function in C 1 (J) such that Ψ (t) = 0 for any t ∈ J and Ψ(0) = 0. Define f (t) = g(Ψ(t)). Then Proof. First, let us observe that Deriving both sides of this relation and dividing by Ψ (t), by using Equation (9), we have Finally, by substituting f (t) − f (0) and g(t) − g(0) in place of f and g and using Equation (10) we conclude the proof.
This proposition leads us to easily give the solution for the relaxation equation whenever Ψ(0) = 0. Indeed, if we define z(t) as the solution of the relaxation equation for the Caputo derivative, hence z(t) = y 0 E ν (at ν ), and y(t) = z(Ψ(t)), we have, by the previous proposition thus y(t) is the solution of Equation (11).

Stochastic Linear fractional-integral Equations with Constant Coefficients and Gaussian Solutions
From now on let us fix a complete filtered space (Ω, F , F t , P).
In this section, we want to study existence and uniqueness of solutions of stochastic linear fractional-integral equations in the form where y 0 , a, b ∈ R, a = 0, and G(t) is a given F t -adapted Gaussian process. From now on, as shorthand notation, let us denote Gaussian process with a.s. r-Hölder continuous paths}.

Remark 1.
Obviously, for any f , g ∈ C 0 (J) and Z ∈ G(J), we have f Z + g ∈ G(J).

The fractional-integral of a Gaussian Process
First, one could ask if the fractional-integral of a Gaussian process is still a Gaussian process. Concerning this problem, we have the following Lemma. Lemma 1. Let Z ∈ G(J) for some time interval J and define Z ν (t) := I ν t Z for t ∈ J. Then Z ν ∈ G(J). Moreover, if J = [0, T] for some T > 0, then Z ν ∈ G ν (J).

Proof. Let us consider
and recall that P(Ω \ A) = 0. Fix ω ∈ A and observe that I ν t Z(·, ω) is well-defined and continuous (in t). We need to show that it is a F t -adapted Gaussian process. Let us define for any ω ∈ A where δ ε n = t−ε n and s j,n = j − 1 2 δ ε n , with j = 1, . . . , n. Hence we have that almost surely Since Z is F t -adapted and with a.s. continuous paths, it is progressively measurable and then, for any n ∈ N and j = 1, . . . , n, Z(s j,n ) is F s j,n -measurable and thus, being s j,n < t, F t -measurable. Hence the variable ∑ n j=1 (t − s j,n ) ν−1 Z(s j,n )δ n is F t -measurable for any n ∈ N and so it is its limit as n → +∞, concluding that for any ε > 0, I ν,ε t Z is F t -measurable. Now let us consider, for k ∈ N, ε k = 1 k . Let us observe that for fixed ω ∈ A we have, for τ ∈ (0, t), which is a L 1 (0, t) function. Hence, we have, by Lebesgue dominated convergence theorem, that . . , t m ) ∈ J m and let us consider the random variable Z ε given by As before, if we define, for fixed i ≤ m and n ∈ N, δ ε i,n := t i −ε n and s i,j,n := j − 1 2 δ ε i,n for j = 1, . . . , n, we have that, by definition of Riemann integral, for any ω ∈ A. Hence we have, for any ω ∈ A, is Gaussian for any n ∈ N. Hence Z ε is almost surely limit of Gaussian random variables, hence it must be Gaussian.
As before, if we consider is almost surely limit of Gaussian random variables and must be Gaussian itself. The arbitrariness of (a 1 , . . . , a m ) ∈ R m and m ∈ N gives us the fact that Z ν is a Gaussian process.
Finally, suppose that J = [0, T] and let us consider t 1 , t 2 ∈ J. Suppose, without loss of generality, that t 1 < t 2 and set concluding the proof.
Let us remark that fractionally integrated Gauss-Markov processes have been also studied in [33].

Compatibility between Fractionally Integrated Gaussian Processes
Now we want to study the behavior of a fractionally integrated Gaussian process with respect to other Gaussian processes. To do this let us first give the following shorthand notation. Definition 7. Let Z(t) and G(t) be two F t -adapted Gaussian processes with a.s. continuous paths. We say that Z(t) and G(t) are compatible if for any n, m ∈ N, any (a 1 , . . . , a n , a n+1 , . . . , a n+m ) ∈ R n+m and any (t 1 , . . . , t n , t n+1 , . . . , t n+m ) ∈ R n+m the random variable is still a Gaussian random variable. This obviously implies that Z(t) + G(t) ∈ G(J). Let us denote, for any G ∈ G(J), It is obvious that if Z(t) and G(t) are independent F t -adapted Gaussian processes with a.s. continuous paths, then Z(t) and G(t) are compatible. Now let us show the following Lemma. ω)) is continuous} and recall that P(Ω \ A) = 0. Fix ω ∈ A and observe that Z ν (t) ∈ G(J) by the previous Lemma and that for ω ∈ A, Z ν (·, ω) is continuous. Thus, we have that t → Z ν (t, ω) + G(t, ω) is continuous for any ω ∈ A. Moreover, since Z ν and G are F t -adapted, Z ν + G is F t -adapted. Now we need to show the compatibility property.

Proof. Let us consider
Let us fix m 1 , m 2 ∈ N, (t 1 , . . . , t m 1 , t m 1 +1 , . . . , t m 2 ) ∈ J m 1 +m 2 and (a 1 , . . . , a m 1 , a m 1 +1 , . . . , a m 2 ) ∈ R m 1 +m 2 and let us define the random variables Let us work on the second one. Fix ω ∈ A. Thus, by recalling the definition of δ ε i,n and s i,j,n for i ≤ m and j = 1, . . . , n given in the previous lemma, we have that hence we have that almost surely where on the RHS we have Gaussian random variables since Z(t) and G(t) are compatible. Hence Z ε is a Gaussian random variable. Moreover, if we define ε k = 1/k for k ∈ N, one has that Z = lim k→+∞ Z ε k almost surely, thus Z is a Gaussian random variable and then Z ν and G are compatible.

Remark 2.
Obviously, for any f , g ∈ C 0 (J) and Z ∈ G(J, G), we have f Z + g ∈ G(J, G). Moreover, for any f , g ∈ C 0 (J) and Z ∈ G(J, G), we also have f Z + gG ∈ G(J, G).

Main Result
Now we are ready to show an existence and uniqueness result in G(J) for the solution of Equation (12) in the fashion of [6] (Theorem 3.3).

Proof. Let us consider
and recall that P(Ω \ A) = 0. Fix ω ∈ A and define the operator A ω : Let us show that A ω is well-posed, i.e., A ω f is continuous. To do this consider t 1 , t 2 ∈ J and suppose, without loss of generality, that t 1 < t 2 . Then, we can set t 2 = t 1 + δ. We have hence, being G(·, ω) continuous, sending δ → 0 + , we have lim δ→0 1 p e γt and then Taking the maximum, we have Thus, we can choose γ > 0 big enough to have With this choice of γ, we have that A ω is a contraction and thus admits a unique fixed point (see [34], Theorem 3.1): let us denote it as Y ν (·, ω).
Moreover, let us consider the sequence (for fixed ω ∈ A) This sequence is such that f n (·, ω) → Y ν (·, ω) in C 0 by contraction theorem (see [34]). Now let us define a stochastic operator A. For ω ∈ A, let us define it as for any stochastic process such that t → f (t, ω) is continuous for any ω ∈ A, while for ω ∈ A let us complete it as we wish, since Ω \ A is a null set.
We can re-interpret our sequence as a sequence of stochastic processes given by Now let us observe that f 0 is a (degenerate) F t -adapted Gaussian process with a.s. continuous paths. For f 1 , we have that which obviously belongs to G(J, G) by Remark 2. Let us suppose that f n−1 ∈ G(J, G). By using Remark 2 and Lemmas 1 and 2, we have that f n ∈ G(J, G). Hence we have that for any n ∈ N, where the limit is in the a.s. sense, thus it is easy to see that Y ν ∈ G(J) (a.s. continuity of the paths follows from the continuity of t → Y ν (t, ω) for ω ∈ A, since Y ν (·, ω) are fixed points of A ω ). Finally, a.s. uniqueness follows from the fact that A ω are contractions for ω ∈ A, hence their fixed point is unique. Now, if G is a.s. ν-Hölder continuous, let us define ω) is ν-Hölder continuous} and let us recall that P(Ω \ A ν ) = 0. Consider ω ∈ A ν and f ∈ C 0 (J) and observe that, from (14), we have In particular, we have that for any f ∈ C 0 (J), A ω f ∈ C ν (J). Almost surely ν-Hölder continuity of the paths of Y ν thus follows from the fact that, for any ω ∈ A ν , Y ν (·, ω) = A ω Y ν (·, ω).

Remark 3.
Let us observe that the fractional-integral operator I ν t is a compact Hilbert-Schmidt operator in L 2 ([0, T]) (see, for instance [35]) if ν > 1 2 . Indeed, the integral kernel k(t, τ) := (t − τ) ν−1 χ (0,t) (τ) (where χ (0,t) is the indicator function of the interval (0, t)) is such that In such case, one can use the structure of the equation to show that there exists a unique Gaussian solution. Setting for instance a = 1 and b = 0, we have for fixed ω ∈ Ω where I is the identity operator; hence we have ω)).
In such a way, for ν > 1 2 , one has the characterization of the solution Y as and then Y is given by a linear operator applied to a Gaussian process, hence it is Gaussian.

Speed of Convergence
We could also investigate the speed of convergence of the sequence f n defined in Equation (15) to Y ν . The following proposition is an easy consequence of the contraction theorem.

Proof.
Fix ω ∈ A (where the set A is defined in Equation (13)). By contraction theorem (see [34]) we have, since L is the Lipschitz constant of A ω : Now let us recall that for any function f ∈ C 0 γ (J) thus, we have Now let us recall that f 0 (t, ω) = 0 and f 1 (t, ω) = y 0 + G(t, ω) = G(t, ω), thus we have .

The Mean of Y ν
Let us introduce another class of Gaussian processes We want to investigate the mean of Y ν , solution of (12), when G ∈ G 0 (J). We have the following result. Proposition 3. Fix G ∈ G 0 (J) and let us suppose that y ν (t) hence y ν (0) = y 0 . Now, let us notice that Hence we can use Fubini's theorem to achieve y ν (t) = y 0 + I ν t (ay ν + b).
Rearranging the equation and applying I 1−ν t on both sides we have Since on the RHS we have a C 1 function, we can differentiate both terms and use (6) and (7) to achieve d ν y ν dt ν (t) = ay ν (t) + b.
The proof of such result is analogous to the previous one.

The Choice of G: Some Examples
In this section, we will give some example concerning the choice of the process G ∈ G(J). Actually, these kinds of equations are noisy versions of the Cauchy problems d ν y ν dt ν = ay ν + b y ν (0) = y 0 (17) and the choice of the noise depends on the choice of G ∈ G(J). Moreover, if we take G ∈ G 0 (J) and E[Y ν (t)] is in L ∞ (J), we are considering a process with an assigned mean and we can modulate covariance by changing G. Let us give some examples.

Brownian Motion and White Noise
A first simple case is given by choosing G(t) as a Brownian motion W(t) on (Ω, F , F t , P). Concerning the regularity of the solution Y ν of Equation (12), we have the following result. (12) ν ∈ (0, 1/2) and G = W a Brownian motion on (Ω, F , F t , P). Then the solution Y ν of Equation (12)  Actually, we could imagine writing (only formally) our integral equation in differential form. We have, by (formally) using the relation (6),

Corollary 1. Consider in Equation
Writing in this way, we can see what the role is of W(t): it works like a white noise introduced in Equation (17).

Fractional Brownian Motion and Fractional White Noise
We could also choose G(t) to be a fractional Brownian motion W H (t) with Hurst index H ∈ (0, 1), introduced in [36]. We will not focus on the features of such process, but for them we refer to [37,38]. Let us recall that the definition of W H already involves fractional integrals. Indeed, for ν ∈ (0, 1), we can define the operators − I ν t and − D ν t as for any suitable function f : R → R. In such case, if H ∈ 1 2 , 1 and we fix ν = H − 1 2 , then there exists a (normalizing) constant C H such that Concerning the regularity of the solution Y ν of Equation (12), we have the following result: Consider H ∈ (0, 1) and let (in Equation (12)) G = W H be a fractional Brownian motion on (Ω, F , F t , P) with Hurst index H; then, if ν ∈ (0, H), the solution Y ν of Equation (12) belongs to G ν (J) ∩ G(J, W H ) for any J = [0, T].
Such corollary is linked to the fact that the paths of W H (t) are γ-Hölder continuous for any γ < H, as shown in [38].
As before, we could formally write Equation (12) in differential form, by using Equation (6), to achieve where the fractional white noise dW H must be carefully interpreted. Thus, we have that our equation is a perturbation of (17) with a fractional white noise. For ν = 1, we obtain the fractional Ornstein-Uhlenbeck process ( [39,40]).

Ornstein-Uhlenbeck Process and Colored Noise
We get another example by choosing G(t) to be an Ornstein-Uhlenbeck process U(t), solution of for some λ, µ ∈ R, σ > 0 and W a Brownian motion on (Ω, F , F t , P). In such case, we have the following regularity result: Corollary 3. Consider, in Equation (12), ν ∈ 0, 1 2 and G := U an Ornstein-Uhlenbeck process on (Ω, F , F t , P). Then, the solution Y ν of Equation (12) belongs to G ν (J) ∩ G(J, U) for any J = [0, T].
As before we can write the differential form of the equation obtaining that, by using Equation (18), becomes thus observing the effect of a colored noise (for ν = 1 see, for instance, [41,42]). Eventually, we could also use a fractional Ornstein-Uhlenbeck U H in place of U, obtaining the following regularity result:

Fractional Itô Integral
There is another particular choice of G that can be done. Let us suppose that ν ∈ 1 2 , 1 and observe that, for fixed t > 0, the function (t − τ) ν−1 is in L 2 (0, t). Thus, the following process is well-defined and belongs to G 0 (J) for any J = [0, T]: where W is a Brownian motion on (Ω, F , F t , P) and the integral must be interpreted in the Itô sense. With this noise, Equation (12) is the integral version of a Caputo-fractional stochastic differential equation, as studied for instance in [24,25]. For such equations, closed form of the solutions can be obtained via a variation of constant formula, as shown in [24]. In this particular case it is known that y ν (t) := E[Y ν (t)] is a continuous function.

Stochastic Models for Fractional Gompertz Curves
In this section, we will construct two classes of stochastic models for fractional Gompertz curves: one for the fractional Gompertz curve given in [21], the other for the improved one given in [22]. Moreover, we will introduce a third model that combines the two previous approaches. However, let us first recall how the classical stochastic Gompertz model is constructed.

The Stochastic Gompertz Model
Here we will recall some basics of the stochastic Gompertz model. We will follow the lines of [43]. Let us consider a stochastic process X(t) solution of the following stochastic differential equation where x 0 > 0 is a constant, W(t) is a Brownian motion (with respect to the filtration F t ) and α, β > 0 are the growth parameters we defined in Section 1.
Now, if we define the process Y(t) = log X(t) x 0 , a simple application of the Itô formula leads to which is the Stochastic Differential Equation of an Ornstein-Uhlenbeck process. It will be useful to write such equation in integral form In particular, Y(t) is a Gaussian process and thus X(t) = x 0 e Y(t) is a log-normal process. Moreover, it is easy to see that, setting y(t) = E[Y(t)], we have dy dt (t) = −βy(t) + α y(0) = 0 that is the second equation of (2). However, since (20) is non-linear, we cannot conclude that the mean of X(t) solves the Gompertz Equation (1). Let us then denote with x(t) the median of X(t), i.e., for each t > 0 Since X(t) is absolutely continuous for t > 0, x(t) is the unique solution of the equation (in z) It is well known that the median of a log-normal variable coincides with the exponential of the mean of its logarithm, i.e., x(t) = x 0 e y(t) .
In particular this shows that x(t) solves Gompertz Equation (1).
For this reason, we can consider X(t) a stochastic version of the Gompertz curve: the deterministic model is recovered via the median of the process, which is an observable function that describes the macroscopic behavior of the process.
Following this line, we are searching for a log-normal (or a Gaussian process) such that the median (or the mean) is a fractional Gompertz curve.

A Stochastic Model for the Fractional Gompertz Curve
Let us first obtain the stochastic model for the fractional Gompertz curve given in [21]. Recalling the construction of the classical stochastic Gompertz model, we are searching for a log-normal stochastic process whose median is a fractional Gompertz curve.
To do this, fix a time horizon T > 0 and a time window J = [0, T]. Let us consider the following stochastic linear fractional-integral equation: (23) can be recognized as a stochastic version of Equation (3). Indeed, the latter can be written in integral form as and then Equation (23) follows by adding a noise G(t). In particular, such equation follows as a generalization of Equation (22) by substituting the classical integral with the fractional one and the white noise with a general Gaussian one.
A natural choice for G(t), if ν ∈ 1 2 , 1 , is the one given in Section 4.4 by Equation (19). In particular, if we set Equation (23) can be formally seen as a fractional version of Equation (21). Now, let us define the process X ν (t) = X 0 e Y ν (t) with Y ν (t) defined in (23). We have the following result.

Proposition 4.
The process X ν (t) is a log-normal process. Moreover, its median x ν (t) is given by i.e., is a fractional Gompertz curve.
Proof. The fact that X ν (t) is a log-normal process follows from the fact that Y ν ∈ G(J). Moreover, since it is a log-normal process, we have where y ν (t) = E[Y ν (t)]. By Proposition 3, we know that y ν (t) is solution of Equation (3) and we conclude the proof.

A Stochastic Model for the Improved Fractional Gompertz Curve
Let us obtain a stochastic model for the improved fractional Gompertz curve. To do this, let us recall that the improved fractional Gompertz curve is solution of In this case, referring to Definition 6, let us choose Ψ as to rewrite Equation (24) as follows: . As before, for ν ∈ 1 2 , 1 , a natural choice for the noise is given by Now let us define X ν (t) := Y ν (Ψ(t)). We have the following result. Proposition 5. X ν (t) is a Gaussian process whose mean x ν (t) := E[X ν (t)] is given by Proof. Recalling that y ν (t) = E[Y ν (t)], we have, from Proposition 3, that y ν is solution of Moreover, recalling also that X ν (t) = Y ν (Ψ(t)), we know that where Ψ(t) = 1 β (1 − e −βt ). Let us observe that x ν (0) = y ν (Ψ(0)) = x 0 . Finally, by Proposition 1, concluding the proof.

A New Fractional Model and Its Stochastic Counterpart
Now we want to give a fractional model that takes into account both the fractional Gompertz curve and the improved fractional Gompertz curve. To do this, we will suppose an a priori form for the growth rate. This has been done for instance in [44], with the introduction of a growth model that takes into account both the Gompertz and the Korf dynamics. In general, let us consider as a starting point a model of the form dx dt (t) = r(t)x(t), for some growth rate function r(t). Equation (25) can be solved and solution is given by In our case, let us consider as growth rate the function for some ν 1 ∈ (0, 1), where E ν,θ (t) is the two-parameters Mittag-Leffler function (see, for instance, [45]) defined as Let us also denote with x ν 1 the solution of Equation (25) where r ν 1 (t) is used in place of r(t). In this case, we can explicitly calculate the integral in Equation (26). Indeed, let us denote R(t) = E ν 1 (−t ν 1 ). We have thus, differentiating the series term by term, we have Now let us consider by using Equation (27) we have Thus, we have that t 0 r ν 1 (t)dt = αΨ ν 1 (t).
By substituting this last integral in Equation (26) for r ν 1 and writing Ψ ν 1 explicitly, we achieve which is actually the fractional Gompertz curve given in [21]. Indeed, recalling that for the fractional Gompertz curve we had x ν 1 (t) = x 0 e y ν 1 (t) where x ν 1 (t) was solution of the equation Since y ν 1 (t) = αΨ ν 1 (t), we have shown that y ν 1 (t) = r ν 1 (t) and we achieve Equation (25). Thus, we can conclude that equation defines the fractional Gompertz curve. Here we already introduced a first degree of fractionalization: now let us introduce another fractional timescale. To do this, let us work as in [22] and let us consider a fractional generalization obtained by introducing the fractional derivative with respect to Ψ ν 1 . So our new fractional Gompertz curve will be defined as the solution x ν (where we denote ν = (ν 1 , ν 2 ) ∈ (0, 1) 2 ) of the fractional Cauchy problem x ν (0) = x 0 that, being a relaxation equation for the Caputo derivative with respect to Ψ, can be explicitly solved as This new fractional Gompertz curve exhibits two degrees of fractionality: one given by the fact that we chose the growth function to be the one of the fractional Gompertz curve, the other from the fact that we introduced a fractional Caputo derivative (with respect to the integral of the growth rate) in the corresponding time in-homogeneous relaxation equation. This also leads to the possibility of considering two different fractional timescales: one for the population, the other for the growth rate. Concerning the stochastic model for such fractional Gompertz curve, fix ν = (ν 1 , ν 2 ) ∈ (0, 1) 2 and let us consider Y ν 2 as the solution of such that y ν 2 (t) = E[Y ν 2 (t)] is in L ∞ (J). As we already stated, if ν 2 ∈ 1 2 , 1 , we could consider Finally, let us define X ν (t) := Y ν 2 (Ψ ν 1 (t)). We have the following result. Proposition 6. X ν (t) is a Gaussian process whose mean x ν (t) := E[X ν (t)] is given by The proof of such proposition is analogous to the one of Proposition 5.

Conclusions
In this paper, we have given some methods to construct stochastic fractional Gompertz models by using stochastic linear fractional equations with Gaussian driving processes. The choice of a Gaussian driving process is linked to the necessity (in the first class of models introduced in Section 5.2) to preserve the lognormality of the Gompertz model. In Sections 5.3 and 5.4 we obtain stochastic fractional Gompertz models with Gaussian one-dimensional law. These were actually only exemplifications. Indeed, one can use the construction method given in Sections 5.3 and 5.4 to obtain Gaussian stochastic models for general growth models of the form x ν (0) = x 0 (28) where Ψ (t) = r(t) α for some growth rate r(t) > 0. One could also try to substitute the operator I ν,Ψ t in place of I ν t in Equation (12). In such a case one could show, by similar arguments, the existence and uniqueness of a Gaussian solution and then use the construction given in Section 5.2 to obtain a log-normal stochastic model for (28).
Concerning possible applications, it has been already observed in [21,22] that fractional Gompertz models are more appropriate than classical ones to describe some phenomena such as tumor growth (concerning the model in Section 5.2), dark fermentation and other fermentation phenomena (concerning the model in Section 5.3). In this paper we provided a method to introduce noise (due to eventual unpredictable variables in the environment) in such a way that a macroscopic observable function still preserves such laws. Concerning the choice of the noise, it depends on the autocorrelation one wants to introduce in the model. For instance, if one wants to introduce a long-range (or short-range) correlated noise, one could use a fractional Brownian motion as driving Gaussian process, while if a delta-correlated noise is needed one could use a classical Brownian motion as driving process.
Finally, we want to recall that our aim was to introduce some construction methods that could lead to log-normal or normal stochastic models for general fractional growth processes (as the ones in Equation (28)) with a general Gaussian noise, in order to provide a wide range of models that could be possibly useful in future applications.
Author Contributions: Conceptualization, G.A. and E.P.; methodology, G.A. and E.P.; writing-original draft preparation, G.A. and E.P.; writing-review and editing, G.A. and E.P. All authors have read and agreed to the published version of the manuscript.
Funding: This research is partially supported by MIUR-PRIN 2017, project "Stochastic Models for Complex Systems", no. 2017JFFHSH.