Generalized Fractional Calculus for Gompertz-Type Models

: This paper focuses on the construction of deterministic and stochastic extensions of the Gompertz curve by means of generalized fractional derivatives induced by complete Bernstein functions. Precisely, we ﬁrst introduce a class of linear stochastic equations involving a generalized fractional integral and we study the properties of its solutions. This is done by proving the existence and uniqueness of Gaussian solutions of such equations via a ﬁxed point argument and then by show-ing that, under suitable conditions, the expected value of the solution solves a generalized fractional linear equation. Regularity of the absolute p -moment functions is proved by using generalized Grönwall inequalities. Deterministic generalized fractional Gompertz curves are introduced by means of Caputo-type generalized fractional derivatives, possibly with respect to other functions. Their stochastic counterparts are then constructed by using the previously considered integral equations to deﬁne a rate process and a generalization of lognormal distributions to ensure that the median of the newly constructed process coincides with the deterministic curve.


Introduction
In the context of population dynamics, the Gompertz curve represents one of the most adaptable models of population growth with variable rate. Precisely, Gompertz proposed the well-known curve to model population growth under the assumption that the mortality rate grows exponentially with the age (see [1]). After this, the Gompertz curve has been shown to be a quite useful model for phenomena that exhibit an intrinsic ageing effect, as for instance tumour incidence (see [2]). Gompertz-type models are not limited to phenomena involving human age. For instance, such kind of models have been widely used in the context of tumour growth (see [3]), as the growth of the radius of a multicell spheroid is influenced by the inhibition effect of its necrotic core (see [4]), which grows together with the tumour itself.
Different generalizations of the Gompertz curve have been presented in literature. An example is given by the generalized logistic curve (see [5]), which covers not only the Gompertz curve, but a wide family of growth curves depending on the choice of the exponents of the logistic equation and the rate function. On the other hand, to include background noise effects, stochastic generalizations of the Gompertz curve, obtained via diffusion processes, have been introduced in [6] and then further extended to nonhomogeneous diffusions in [7], while first passage time problems for them have been studied, e.g., in [8,9]. Such models have been widely used, for instance, to study the effects of therapy on tumours (see [10] and references therein). The importance of the stochastic interpretation of the Gompertz curve is underlined in [11], where the crucial role of the noise in growth phenomena is highlighted.
Among the generalizations of Gompertz-type models, we also find fractional order Gompertz curves. Fractional calculus has been applied to a lot of different fields (see [12] growth phenomena is also mandatory to describe possible random fluctuations, as done in [6,7] for the classical Gompertz curve and in [29] for the fractional Gompertz curves.
The problem we address in this paper is twofold: on one hand, we want to extend the definition of fractional Gompert curves to cover a wider range of memory kernels; on the other hand, we want to introduce noise in the aforementioned models in a coherent way. Thus, here we consider a further generalization of the construction presented in [29] to the case of fractional orders expressed in terms of complete Bernstein functions. In particular, we first construct some deterministic generalized fractional Gompertz curves and then we introduce their stochastic counterparts. Differently from [29], we consider a generalization of lognormal processes to guarantee that the stochastic Gompertz curves remain positive and, at the same time, their medians represent the deterministic ones. The paper is structured as follows: • In Section 2 we give some preliminaries on Bernstein functions; • In Section 3 we introduce generalized fractional integrals and derivatives. In particular we give some chain rules involving generalized fractional derivatives with respect to other functions and we study the properties of the eigenfunctions of the defined operators; • In Section 4, we study the properties of the solutions of some linear stochastic equations involving generalized fractional integrals. In particular, we show that such equations admit a unique Gaussian solution and that, under some suitable assumptions on the noise process, its expectation solves a linear generalized fractional differential equation; • In Section 5, we use the operators introduced in Section 3 to construct generalized fractional Gompertz curves and the Gaussian processes given in Section 4 to determine their stochastic counterparts; • In Section 6 we present a short summary of the main results of the paper and we compare them with some pre-existing literature. Finally, some highlights for future works are given.

Preliminaries and Notation
In this section we provide some preliminary definitions and lemmas.
The convex cone (see [34] [Corollary 3.7]) of Bernstein functions will be denoted as BF . A Bernstein function Φ is said to be special if and only if the conjugate function is still a Bernstein function. The class of special Bernstein functions will be denoted as SBF .
To characterize Bernstein functions, let us recall the following theorem, known as Lévy-Khinchine representation theorem (see [34] [Theorem 3.2]). Theorem 1. A function Φ : R + → R belongs to BF if and only if there exist two constants a Φ , b Φ ≥ 0 and a non-negative measure ν Φ on R + such that and The measure ν Φ is called the Lévy measure of Φ. Vice versa, any triplet (a Φ , b Φ , ν Φ ), where a Φ , b Φ ≥ 0 and ν Φ is a non-negative measure on R + satisfying condition (1), defines a unique Φ ∈ BF via Equation (2).
The constants a Φ and b Φ are called respectively the killing and the drift coefficient of Φ. We will denote ν Φ (t) = ν Φ (t, +∞). However, we need Bernstein functions whose Lévy measure is more regular. To do this, we refer to the following definition.

Definition 2.
A Bernstein function Φ is said to be complete if its Lévy measure ν Φ (dt) admits a completely monotone density, that is, We denote the convex cone (see [34] [Corollary 7.6]) of complete Bernstein functions as CBF .
Bernstein functions can be recognized as Laplace exponents of particular Lévy processes. Indeed, let us recall the following definition.

Definition 3.
A subordinator σ(t) (see [43] [Chapter I I I]) is a non-decreasing Lévy process. Given a subordinator σ(t) and a positive constant a > 0, the process where τ a is an exponentially distributed random variable, with rate a, independent of σ(t), is called a subordinator killed at rate a.
Concerning the link between subordinators and Bernstein functions, we have the following Theorem (see [34] [Theorem 5.2]).
Since we focus on Bernstein functions, we will usually denote a subordinator by σ Φ (t), referring to the fact that its Laplace exponent is given by Φ.
For any subordinator, the following occupation measure can be defined.
is the Borel σ-algebra of R + 0 and 1 A is the indicator function of the Borel set A. We will denote the distribution function of the potential measure U Φ (t) := U Φ (0, t) and we will usually refer to it directly as a potential measure.
Moreover, we can define a right-continuous inverse for the process σ Φ (t).

Remark 1. As shown in
Once we have defined the inverse subordinator, by using the fact that σ Φ and L Φ are increasing, we have Concerning special Bernstein functions, the associated potential measure is almost absolutely continuous, except at most for a jump in 0, as stated in the following theorem (see [34] [Theorem 10.3]). Theorem 3. Let Φ ∈ BF . Then Φ ∈ SBF if and only if there exists a non-negative and non-increasing function u Φ (t) such that In particular, if b Φ > 0 or b Φ = 0 and ν Φ (R + ) = +∞, then U Φ is absolutely continuous with density given by u Φ , called the potential density of Φ.

Remark 2.
Let us emphasize that, for fixed Φ ∈ SBF , denoting by a Φ , b Φ , ν Φ respectively the killing coefficient, the drift coefficient and the Lévy measure of Φ , it holds (see [34] [Equation (10.9)]): Let us recall that the following inclusion holds: We will work with a specific subset of complete Bernstein functions. Hence, let us introduce the following notation: where ι(λ) = λ is the identity map. Hereafter, we consider the following Assumption. Assumption 1. Φ ∈ CBF and there exist t 0 > 0, C > 0 and β ∈ (0, 1) such that

Remark 3.
The previous Assumption guarantees that, for any T > 0, there exists a constant C such that

Generalized Fractional Integrals and Derivatives
Let us first introduce some generalized fractional integrals.
Definition 6. Set a Banach space (X, | · |). For any T > 0, given a function f : [0, T] → X, we say that f ∈ L 1 ([0, T]; X) if f is Bochner-integrable, that, by Bochner's Theorem (see [44] [Theorem 1.1.4]) is equivalent to asking that f is measurable and . Given a function f : For any function f ∈ L 1 loc (R + 0 ; X) we define the generalized fractional integral of the first kind of f induced by Φ ∈ CBF as where u Φ is the potential density of Φ and the integral is a Bochner integral. We define the generalized fractional integral of the second kind of f induced by Φ ∈ CBF as where ν Φ is the tail of the Lévy measure of Φ and the integral is a Bochner integral. Let us also underscore that, by Remark 2, as Φ ∈ CBF , the operators I Φ and I Φ coincide if and only if (0,+∞) tν Φ (dt) = +∞ or Φ = ι. In general, it holds, for any f ∈ L 1 loc (R + 0 ; X), With the help of the previously introduced operators, we can define the following generalizations of both Riemann-Liouville and Caputo derivatives, introduced first in [32] for complete Bernstein functions and then in [33] in the general case. Definition 7. Set a Banach space (X, | · |). For any function f ∈ L 1 (R + 0 ; X) we define the generalized Riemann-Liouville derivative induced by Φ ∈ CBF , with Φ = ι, as where the integral is a Bochner integral, provided that the involved quantities exist. Given T > 0, we say that a function f : [0, T] → X is absolutely continuous if there exists a function f : (0, T) → X such that f ∈ L 1 ((0, T); X) and and we denote it by f ∈ AC([0, T]; X). We denote by AC loc (R + 0 ; X) the set of functions f : R + 0 → X such that for any T > 0 it holds f ∈ AC([0, T]; X). When X = R, we do not specify X.
For any f ∈ AC loc (R + 0 ; X) we define the generalized Caputo derivative of f induced by Φ ∈ CBF , with Φ = ι, as where the integral is a Bochner integral.

Remark 5.
Let us remark that, if f ∈ AC loc (R + 0 ; X), then one can show that I Φ f ∈ AC loc (R + 0 ; X). Thus, the quantity defined in Equation (7) is well-defined for any f ∈ AC loc (R + 0 ; X). However, if Φ ∈ CBF \{ι}, there exist some functions f ∈ AC loc (R + 0 ; X) such that I Φ f ∈ AC loc (R + 0 ; X). For instance, if X = R, this is the case of f (t) = ν Φ (t). Indeed, by Equation (6) and [34] [Theorem 10.9] we get Such operators are generalizations of the well-known Caputo and Riemann-Liouville fractional derivatives, which are achieved in the case Φ(λ) = λ α with α ∈ (0, 1). Indeed, in this case, In general, for Φ ∈ BF such that a Φ = 0 the operators are defined as Let us consider any f ∈ AC loc (R + 0 ; X). By a simple application of Fubini's theorem we get where we also used [45] [Chapter 6, Theorem 11]. Differentiating (almost everywhere) in both sides of the previous relation we get , for almost any t > 0. (10) With this relation in mind we can extend the definition of generalized Caputo derivative to a (possibly) larger class of functions. Definition 8. Set a Banach space (X, | · |). For any function f : R 0 + → X we define the regularized generalized Caputo derivative induced by Φ ∈ CBF as whenever the right-hand side is well-defined.
Remark 7. Equation (10) justifies the fact that we are using the same symbol of the generalized Caputo derivative.
Concerning the inversion of such operators, let us observe that, as shown in [36] thus we can see the operator I Φ t as the inverse of the generalized Caputo derivative d Φ dt Φ . In the case Φ(λ) = λ α , it holds I Φ = I Φ = I 1−α , where I 1−α is the fractional integral of order 1 − α (see [13] [Chapter 2]), D Φ t = D α t is the Riemann-Liouville fractional derivative of order α and d Φ dt Φ = d α dt α is the Caputo fractional derivative of order α. Here, we also need other generalized fractional operators, that is, the integral and the derivative of a function with respect to an increasing function. The definitions we give are analogous to the ones given in [27] for the Riemann-Liouville fractional derivative and [26] for the Caputo one. Definition 9. Set a Banach space (X, | · |) and consider a strictly increasing function Ψ ∈ C 1 (R + ) ∩ C 0 (R + 0 ). For any measurable function f : R + 0 → X, we define the generalized fractional integral of the second kind of f induced by Φ ∈ CBF with respect to the function Ψ as where ν Φ is the tail of the Lévy measure of Φ and the integral is a Bochner integral, provided that the involved quantities exist. Suppose now Ψ (t) = 0 some t > 0. For any measurable function f : R + 0 → X we define the generalized Riemann-Liouville derivative induced by Φ ∈ CBF , Φ = ι, with respect to Ψ as where the integral is a Bochner integral, provided that the involved quantities exist.
Moreover, for any f ∈ AC loc (R + 0 ; X), we define the generalized Caputo derivative of f induced by Φ ∈ CBF , Φ = ι, as where the integral is a Bochner integral, provided that the involved quantities exist.

Remark 8.
Let us underline that we do not really need Ψ (t) = 0 for all t > 0, but only on the points in which we want to define D Φ,Ψ t . Moreover, we can formally define the Caputo derivative even if Ψ (t) = 0, since Ψ does not play any role in the first equality of Formula (13).
Let us stress that if Φ(λ) = λ α the operators coincide with the ones introduced in [26,27]. Now we want to prove a chain rule, analogous to the one given in [29] [Proposition 1].
It holds 2. If , provided one of the involved quantities exists; 3.
If g ∈ AC loc (R + 0 ; X) and Ψ is locally Lipschitz in R + 0 , then f ∈ AC loc (R + 0 ; X) and it holds Proof. Let us argue for Φ = ι, since the case Φ = ι is trivial.
By the definitions of I Φ,Ψ in Equation (12) and I Φ in Equation (5) we have where we used the change of variables s = Ψ(τ).
Concerning claim (2), it follows from (1) by considering Φ in place of Φ, differentiating both sides of (14) and multiplying by 1 Ψ (t) . Let us prove claim (3). First of all f ∈ AC loc (R + 0 ; X) since it is composition of an absolutely continuous function with a locally Lipschitz one. In particular, it holds , for almost any t > 0 and then where we again used the change of variables s = Ψ(τ). Finally, concerning claim (4), we have, by claims (2) and (3) and the fact that Ψ(0) = 0 (and thus f (0) = g(0)), Remark 9. Let us observe that claims (3) and (4) also hold if Ψ is not locally Lipschitz, but g ∈ C 1 (R + ) ∩ C 0 (R + 0 ). Indeed, both claims directly follow from the fact that in such case Moreover, the last Proposition tells us that the quantity (12) is well defined for any measurable function f : . This is the case, for instance, of f ∈ C 0 (R + 0 ), since, by Continuous Inverse Theorem [46] [Theorem 5.6.5], f • Ψ −1 ∈ C 0 (R + 0 ). Finally, the quantity in Equation (13) is well defined if Ψ is locally Lipschitz in R + 0 whenever f ∈ AC loc (R + 0 ; X) with the property that there exists g ∈ AC loc (R + 0 ; X) such that f (t) = g(Ψ(t)) (this is, for instance, the case in which f ∈ AC loc (R + 0 ; X) and Ψ is bi-Lipschitz). If Ψ is not locally Lipschitz, then the quantity in Equation (13) is still well-defined if the function g defined as above belongs to

Eigenfunctions of the Generalized Fractional Derivatives of Caputo Type
In the following we need to characterize the eigenfunctions of the generalized fractional derivatives of Caputo type that have been previously introduced. First of all, let us recall that the function, is well-defined for any λ ∈ R (see [38] [Lemma 4.1]). Let us also recall the following Proposition (see [38] [Proposition 4.3]).

Remark 11.
As a consequence of [41] In the following we need to extend the definition of e Φ for fixed λ ∈ R to negative values of t. To do this we set that is a continuous monotone function.
e ± Φ (·; λ) can be recognized as the eigenfunction of a particular non-local operator. Indeed, let us define the following operator. Definition 10. Set a Banach space (X, | · |). We say that a measurable function f : , if and only if the function g : R + 0 → X, defined as g(t) = f (−t) for any t ≥ 0, belongs to L 1 loc (R + 0 ; X). For any function f ∈ L 1 loc (R − 0 ; X) we define the right generalized Riemann-Liouville deriva- We say that f : Since for any f ∈ AC loc (R − 0 ; X) it holds we can extend the definition of right generalized Caputo derivative to non absolutely continuous functions via Equation (15), supposed that the function admits a generalized Riemann-Liouville derivative.
. For any function f ∈ L 1 loc (R; X) we define the bilateral generalized Riemann-Liouville derivative induced by Φ ∈ CBF as and the bilateral generalized Caputo derivative induced by Φ ∈ CBF as provided the involved quantities exist.

Remark 13.
The definition of right derivative is given by taking in consideration the integration by parts formula (see, for instance [28]). However, with d Φ,± dt Φ we want to mimic the behaviour of the derivative on the whole real line, thus we need to introduce another − sign on the right derivative.
Obviously, the bilateral derivatives are well-defined on AC ± loc (R; X).

Proposition 3.
Fix Φ ∈ CBF and let g ∈ L 1 loc (R + 0 ; X) and f ∈ L 1 (R − 0 ; X) such that f (t) = g(−t) for any t ≤ 0. The following properties are true: , provided one of the involved quantities exists;

2.
It holds , provided one of the involved quantities exists.
Proof. Let us first observe that where we used the change of variables s = −τ. Differentiating on both sides we get By using (15), we conclude the proof.
As a direct consequence of Propositions 1-3, we obtain the following result.
The following statements hold: We already know, by definition, that sign(λ) e ± Φ (t; λ) is nondecreasing. We want to prove that sign(λ) e ± Φ (t; λ) is strictly increasing. We actually have a stronger result.

Gaussian Solutions for a Linear Stochastic Integral Equation with Constant Coefficients
From now on, let us fix a complete filtered probability space (Ω, F , {F t } t≥0 , P). For a fixed Φ ∈ CBF , we want to exhibit the solution of the following stochastic integral equation where G is a suitable Gaussian stochastic process, Y 0 is a Gaussian random variable independent of G and a, b ∈ R. Before proving an existence and uniqueness result, we need to set some notations.

Properties of Generalized Fractionally-Integrated Gaussian Processes
For any 0 < T ≤ +∞, let us denote by J = [0, T] a time interval with horizon T. If T = +∞, we set J = R + 0 . For T = +∞, we can define the Banach space (C 0 (J), · 0 ) of continuous functions equipped with the supremum norm Let us consider the following class of stochastic processes: and its subclass of Gaussian processes: To study Hölder-regularity properties of the sample paths of the solution of (21), we consider the following subclass of G(J): As a preliminary result, let us show that if we fix 0 < T ≤ +∞ and we integrate a process Z ∈ G(J), we obtain a process in G(J). Lemma 1. Let 0 < T ≤ +∞, Φ ∈ CBF and Z ∈ G(J). Then there exists a set A ⊆ Ω such that P(A) = 1 and the process Z Φ (ω, t) := I Φ Z(ω, ·)(t) is well defined for any ω ∈ A. Moreover, if Assumption 1 is satisfied, then Z Φ ∈ G β (J).

Proof. Let
A = {ω ∈ Ω : Z(ω, ·) ∈ C 0 (J)} and observe that, being Z ∈ G(J), P(A) = 1. Hence, for any ω ∈ A, Z Φ (ω, ·) is welldefined. We omit the proof of the fact that Z Φ ∈ G(J) as it is identical to the first part of the proof of [29] [Lemma 1]. Let us show that Z Φ ∈ G β (J). To do this, fix any T ∈ J, ω ∈ A and define J = [0, T ]. We want to show that Z Φ (ω, ·) is β-Hölder continuous in J . To do this, fix t ∈ J and h ∈ R such that t + h ∈ J . Let us first assume h > 0. We have where we could take the supremum norm since, being ω ∈ A, Z(ω, ·) ∈ C 0 (J ), and we also used the fact that If h < 0 the argument is analogous except for the fact that we control U Φ (t) − U Φ (t + h) ≤ U Φ (|h|), by subadditivity of the potential measure. Thus, in general, we get concluding the proof. Now let us recall the definition of compatibility, as given in [29]. Definition 11. Let G, Z ∈ G(J). We say that Z is compatible with G if the coupled process (Z, G) is a F t -adapted Gaussian process with a.s. continuous sample paths. Let us recall that this implies that Z + G ∈ G(J).
We denote G(J, G) = {Z ∈ G(J) : Z is compatible with G on J}.
Let us give the following Lemma.
We omit the proof, since it is identical to the one of [29] [Lemma 2].

Remark 16.
Let us also stress that any continuous function f : J → R belongs to G(J, G) (considering it as a degenerate stochastic process). Moreover, if Z ∈ G(J) is independent of G, then Z ∈ G(J, G). Finally, if Z ∈ G(J, G) and f i : J → R are continuous functions for i = 1, 2, 3, then

Existence and Uniqueness of a Gaussian Solution
Now we are ready to prove the following Theorem.

Theorem 4.
Consider Φ ∈ CBF satisfying Assumption 1, G ∈ G(R + 0 ) and Y 0 a Gaussian random variable independent of G. Then Equation (21) admits a Gaussian solution Y ∈ G(R + 0 ), in the sense that there exists a Gaussian process Y ∈ G(R + 0 ) and a set A ⊆ Ω with P(A) = 1 such that Moreover, the solution is unique, in the sense that, if Y is another solution of (21), Finally, if G ∈ G β (R + 0 ) for some β ∈ (0, 1], then Y ∈ G min{ β,β} (R + 0 ).

Proof.
First of all, fix T > 0 and define J = [0, T]. On the space C 0 (J) of continuous functions on J define the norm thus (C 0 (J), · γ ) is a Banach space. Let us consider A = {ω ∈ Ω : G(ω, ·) is continuous} (24) and fix ω ∈ A. Define the operator A (T) Let us first show that A (T) ω is well-defined. Consider f ∈ C 0 (J), t ∈ J and δ ∈ R such that t + δ ∈ J. If δ > 0, we have where we used the fact that U Φ (t) − U Φ (t + δ) ≤ 0 and Remark 3. If δ < 0, arguing in the same way, we get where we used, in this case, the fact that U Φ (t) − U Φ (t + δ) ≤ U Φ (|δ|). Being ω ∈ A, this is enough to guarantee that A (T) ω f ∈ C 0 (J).

Now let us show that A (T)
ω is a contraction. To do this, consider f i ∈ C 0 (J), i = 1, 2, and observe that where we used Assumption 1. Now let us observe that 0 > β − 1 > −1, thus we can choose p > 1 such that 0 > p(β − 1) > −1. Let q > 1 be such that 1 p + 1 q = 1. By Hölder's inequality we get Multiplying both sides of the previous inequality by e −γt we get Now let us observe that hence we can choose γ to be big enough to have ω is a contraction on (C 0 (J), · γ ). By the contraction theorem (see [48] [Theorem 3.1]), we know that A (T) ω admits a unique fixed point in C 0 (J) that we denote Y (T) (ω). Now we need to extend (in some sense) this solution to the space C 0 (R + 0 ). For any ω ∈ A, let us define the operator A ω : C 0 (R + 0 ) → C 0 (R + 0 ) as: Let us first show that A ω is well defined. Indeed, if T 1 > T 2 > 0 and t < T 2 we have Now let us show that A ω admits a fixed point. Let us consider the function Let us first show that it is well-posed. Fix T 1 > T 2 > 0 and t ≤ T 2 . Then we have Now let us show that Y(ω) is a fixed point for A ω . To do this, consider any t ≥ 0 and T > t. Then we have Being t ≥ 0 arbitrary we conclude that A ω Y(ω) = Y(ω). Now let us show that Y(ω) is the unique fixed point of A ω . To do this, suppose A ω admits another fixed point Y (ω). Consider any T > 0. Then we have Being T > 0 arbitrary we get Y (ω, t) = Y(ω, t) for any t ≥ 0. Now, for any stochastic process f ∈ S(R + 0 ), let us define the set Then we can define the operator A : where we recall that A is defined in Equation (24). Let us also define Y(ω, ·) ≡ 0 as ω ∈ Ω \ A, so that Y ∈ S(R + 0 ) and A Y = Ω. Now let us show that Y is a fixed point for A.
thus Y is solution of Equation (21). Now let us show that Y is a Gaussian process. By definition, it is sufficient to show that Y (T) are Gaussian processes for any T > 0. Moreover, since P(A) = 1, we can restrict the probability space to (A, F , {F t }, P) without loss of generality. Fix T > 0, set J = [0, T] and consider the sequence of stochastic processes: ω is a contraction for fixed ω ∈ A, we have that f n → Y (T) pointwise with respect to ω ∈ Ω. Moreover f 0 is degenerate, hence f 0 ∈ G(J, G). Let us suppose f n ∈ G(J, G). Then, by Lemma 2 and Remark 16 we have that f n+1 ∈ G(J, G). By induction, we get that f n ∈ G(J, G) for any n ≥ 0 and in particular, being the limit of Gaussian processes, Y (T) ∈ G(J). Since T > 0 is arbitrary we conclude that Y ∈ G(R + 0 ). Now let us show the uniqueness of the solution. Let Y ∈ G(R + 0 ) be another solution of Equation (21) and let A be the set on which (22) holds for Y . Let A = A ∩ A and observe that P(A ) = 1. Let ω ∈ A and observe that that is, the function Y (ω, ·) is a fixed point of A ω . However, we know that A ω admits a unique fixed point Y(ω, ·), thus Y(ω) = Y (ω) for ω ∈ A . In conclusion, we get Now let us show the last part of the statement of the Theorem. Suppose G ∈ G β (R + 0 ) for some β ∈ (0, 1]. Then there exists a set B ⊆ A with P(B) = 1 such that for any ω ∈ B and any T > 0 it holds: for some constant H(ω, T) > 0. In particular, if β < β, Equation (26) Using Y in place of f in the previous inequality and the fact that A thus obtaining, in this case, Y ∈ G β (R + 0 ). This concludes the proof.

Remark 17.
The initial guess f 0 (ω, t) ≡ 0 in Equation (29) can be substituted with any process f 0 ∈ S(J). This tells us, in a certain sense, that the Picard iteration method proposed to solve (21) in J = [0, T] for some T ∈ R + is stable with domain of attraction S(J).

Speed of Convergence of the Iteration Method
We can use again the contraction theorem to estimate the speed of convergence of the iteration procedure we presented in Theorem 4 in the finite-horizon case. Proposition 6. Consider Φ ∈ CBF satisfying Assumption 1, T ∈ R + , J = [0, T], G ∈ G(J) and Y 0 a Gaussian random variable independent of G. Consider A as in Equation (24) and define the operator A : S(J) → S(J) as ω is defined in Equation (25) and A f is defined in Equation (28). Fix f 0 ∈ S(J) and define f n = A f n−1 for n ≥ 1. Consider any p ∈ 0, 1 1−β , q > 1 such that 1 p + 1 q = 1 and γ > 0 big enough to have where C and β are defined in Assumption 1. Let Y ∈ G(J) be the unique Gaussian solution of Equation (21) and define Then, almost surely. As a consequence, for any ε > 0, it holds Proof. First of all, observe that P(A ∩ A f 0 ) = 1 and that A f n = A f 0 for any n ≥ 0. Fix ω ∈ A ∩ A f 0 . We have shown in Theorem 4 (precisely Equation (27)) that, with our choice of p and γ, A ω is a contraction on (C 0 (J), · γ ) with Lipschitz constant L < 1 defined in Equation (30). Then, by the contraction theorem (see [48] [Theorem 3.1]), we know that By using the equivalence relation given in Equation (23) and recognizing G = f 1 − f 0 , we get concluding the proof.
The right-hand side of Equation (32) could be difficult to evaluate explicitly. However, in a particular case, we can provide a less sharp but explicit bound. Corollary 2. With the notation of Proposition 6, suppose f 0 ∈ C 0 (J), Y 0 ∈ R is deterministic and G is a martingale with E[|G(t)| 2 ] = σ 2 (t). Set Then, for any ε > 0 and n > log((1−L)ε)−log(Me γT ) log(L) , it holds Proof. Fix ω ∈ A (recall that A f 0 = Ω) and observe that, by the triangular inequality, Thus, by (31) and recalling the definition of M in Equation (33), we have that almost surely. Hence, for any ε > 0, it holds , so that Since G is a martingale, |G| is a non-negative integrable submartingale. Thus, by Doob's L 2 inequality (see [49] [Chapter II, Theorem 1.7]), concluding the proof.

Estimates on the Moments of Y
Now we want to give some estimates concerning the integrability of the absolute moments of the solution of (21).
Let Y ∈ G(R + 0 ) be the unique solution of (21), as provided in Theorem 4. It being a Gaussian process, we already know that E[|Y(t)| p ] is finite for any p ≥ 1 and any fixed t ∈ R + 0 . First, we want to determine some sufficient conditions under which the function M 1 (t) = E[|Y(t)|] is locally integrable or bounded. Indeed, under some regularity assumptions on the Gaussian noise process G, we have the following result.

Proposition 7.
Consider Φ ∈ CBF satisfying Assumption 1 and G ∈ G(R + 0 ) such that E[|G(·)|] ∈ L 1 loc (R + 0 ). Then, for any T > 0, there exists a constant C > 0 such that for t ∈ [0, T], where E β is the Mittag-Leffler function defined in Equation (20). In particular, if Proof. Let us recall that (21) holds almost surely. In particular, applying the absolute value and the expectation operator on both sides of (21), we have where we also used the triangular inequality.
Fix any T > 0 and observe that f 1 ∈ L 1 (0, T). Then, by the generalized Grönwall inequality in [38] [Theorem 5.1], we get for some constant C > 0, concluding the first part of the proof.
where the right-hand side is an increasing function. Thus, by [38] [Theorem 5.1] we get As T > 0 is arbitrary, we conclude that M 1 ∈ L ∞ loc (R + 0 ).
As a direct consequence of the previous result we have the following Corollary.

Corollary 3.
Consider Φ ∈ CBF satisfying Assumption 1, G ∈ G(R + 0 ) and define that is to say Proof. Let us first stress out that |G(t)| is a folded Gaussian variable for each t ≥ 0, thus and, in particular, E[|G(·)|] ∈ L ∞ loc (R + 0 ). By Proposition 7 we have that M 1 ∈ L ∞ loc (R + 0 ) and then Taking the expectation on both sides of (21) and using Fubini's theorem, justified by inequality (36), we get Observing that m(0) = E[Y 0 ] and taking the Caputo-type derivative on both sides we see that the function m(t) must be a solution of the Cauchy problem (34). Now let us observe that the function is also solution of the Cauchy problem (34), thus, by uniqueness of global solutions (see [38] [Corollary 3.7]) we conclude the proof.
The previous Corollary will play a major role in defining generalized fractional stochastic Gompertz curves.
The arguments we adopted to show that M 1 is locally bounded can be generalized to any absolute p-moment. Indeed, let us consider p ≥ 1 and M p (t) := E[|Y(t)| p ]. We have the following result.
In this particular case, x n 1 F 1 −n, 1 2 , − 1 2 y 2 x is a polynomial in the real variables x, y. Precisely, we have (see [50] [Formula 8.953]) x n 1 F 1 −n, where H 2n is the Hermite polynomial of degree 2n. Let us recall that the Hermite polynomials are defined as H n (x) = (−1) n e x 2 d n dx n e −x 2 , n ∈ N .
Thus, being M 2n (t) a polynomial function of M 2 (t) and m(t), that are both in L ∞ loc (J), we conclude that M 2n (·) ∈ L ∞ loc (J). Finally, if p = 2n, then there exists n ∈ N such that p < 2n and we conclude the proof by Hölder's inequality.

Generalized Fractional Gompertz Curves
Let us first introduce the deterministic generalized fractional Gompertz curves. To do this, let us recall that, in the classical setting, a Gompertz curve x(t) is defined as the unique solution of the non-linear ordinary differential equation, where a, b, x 0 > 0, that is to say Its rate function is defined as and it is the unique solution of With this auxiliary function, we can recast Gompertz Equation (38) as the following system of bilinear ordinary differential equations Let us highlight that y (t) = ae −bt > 0 for any t > 0, thus we can recast again (39) to achieve that is, we have distinguished the two main components of a Gompertz curve by two linear equations: • The rate function y(t) satisfies a non-homogeneous linear differential equation with initial condition y(0) = 0; • The Gompertz curve itself x(t) satisfies a linear equation with respect to the operator 1 y (t) d dt ; precisely, it is an eigenfunction with eigenvalue 1 of such operator. Now that we have two linear equations in system (40), one can apply a fractionalization procedure on each of these equations. In [24] the authors use a Caputo fractional derivative in place of the standard one in the equation of the rate function, obtaining the system, where α ∈ (0, 1). In this case, the rate function and the curve are given by On the other hand, in [25] the authors propose a Gompertz model in which the fractionalization procedure is applied to the curve equation, obtaining where α ∈ (0, 1). In this case the rate function and the curve are given by Let us recall that the function is strictly increasing by Corollary 1 (since it coincides with e ± Φ (t; 1) as Φ(λ) = λ α , see [40]) and then we can define the function log α as the inverse of e ± α . Thus, the relation that links the rate function with the curve is given by coherently with the chain rule formula given in Proposition 1.
In [29], both approaches are used, obtaining the system where (α 1 , α 2 ) ∈ (0, 1] 2 , and then the rate function and the curve are given by Let us emphasize that if α 1 = 1 we obtain (41) and if α 2 = 1 we obtain (42). Here we want to extend this approach to the case of any complete Bernstein function.
Let us consider (Φ 1 , Φ 2 ) ∈ CBF 2 := CBF × CBF . We define the generalized fractional Gompertz system of generalized fractional orders (Φ 1 , Φ 2 ) as the following system of equations We obtain an explicit expression of both the rate function and the curve by using Propositions 2 and 4 respectively, and we also have the relation In particular, if we choose Φ i (λ) = λ α i , with α i ∈ (0, 1], i = 1, 2, we obtain again Equation (43).

Generalized Fractional Stochastic Gompertz Curves
Now that we have defined a generalization of the (deterministic) fractional Gompertz curves introduced in [24,25,29], we want to construct the respective stochastic versions. Here, we will follow a slightly different route compared to [29].
Indeed, let us first give a generalization of the class of lognormal processes, by using the functions log Φ defined before. Definition 13. Let {X(t), t ≥ 0} be a stochastic process with X(t) > 0 almost surely for any t ≥ 0 and let Φ ∈ CBF . We say that X is a log Φ -normal process if the process Y(t) = log Φ is a Gaussian process.

Remark 19.
If Φ = ι, the definition coincides with that of the lognormal process. Now we are ready to define the stochastic version of the generalized fractional Gompertz curve.
and fix x 0 > 0. We call generalized fractional stochastic Gompertz curve, with generalized fractional orders (Φ 1 , Φ 2 ) and noise process G, the process and we call Y(t) its rate process.

Remark 20.
It is easy to see that, by definition, X(t) is a log Φ 1 -normal process.
Let us first show how such process generalizes the stochastic Gompertz process. Usually, a stochastic Gompertz process is defined via the non-linear stochastic differential equation where W(t) is a standard Brownian motion and k = x 0 e a/b is called carrying capacity. It can be shown by a direct application of Itô's formula (see [11]) that the process x 0 is an Ornstein-Uhlenbeck process, solution of the stochastic differential equation Thus we can equivalently define a stochastic Gompertz process as the process that is the integral formulation of (47). With this equivalent definition, we can recognize the (classical) stochastic Gompertz curve as a particular case of Definition 14, with Φ 1 = Φ 2 = ι and G = √ 2aW. Moreover, by choosing Φ 1 = ι and Φ 2 (λ) = λ α 2 , for some α 2 ∈ (0, 1], in Definition 14, we obtain the fractional stochastic Gompertz curve defined in [29] [Section 5.2]. On the other hand, the curves defined in [29] [Sections 5.3 and 5.4] are Gaussian processes and then they cannot be constructed by means of Definition 14.
A first simple but crucial property of the generalized fractional stochastic Gompertz curve is shown in the following Proposition.
implying y = m(t), that is a contradiction. Hence (48) holds. Now we can show the link between the stochastic growth curves given in Definition 14 and the deterministic ones defined in the previous Subsection.
Theorem 5. Let X(t) be a generalized fractional stochastic Gompertz curve with generalized fractional orders (Φ 1 , Φ 2 ) ∈ CBF 2 , noise process G ∈ G(R + 0 ) and starting point x 0 > 0 and denote by Y(t) its rate process. Let M(t) be the median of X(t) with M(0) = x 0 and m(t) = E[Y(t)]. Then M(t) is a generalized fractional Gompertz curve of generalized fractional orders (Φ 1 , Φ 2 ) with rate function m(t), that is, (M(t), m(t)) is the unique solution of Equation (44).
Proof. Let us first recall that E[G(t)] ≡ 0 and E[G 2 (t)] ∈ C 0 (J) by Definition 14. Hence, by Corollary 3, we know that m(t) is the unique solution of that is to say Now let us stress out that, being Φ 2 ∈ CBF , the function t → e Φ 2 (t; −b) is completely monotone (as shown in [36] [Theorem 2.1]) thus, in particular, m(t) is strictly increasing and differentiable for all t > 0, with m (t) > 0.
By Proposition 3 it holds and then, by Proposition 4, we know that M(t) is the unique solution of concluding the proof.

Conclusions
In this paper, we used the theory of complete Bernstein functions and the tools from generalized fractional calculus to extend the fractional Gompertz curves introduced in [24,25,29]. Precisely, the classical Gompertz equation is decomposed in the rate equation and the curve equation and then a fractionalization (or, in this case, nonlocalization) procedure is applied to both, obtaining the generalized fractional Gompertz curves. For their stochastic counterparts, we first studied in Section 4 a linear integral equation that plays the role of the equation of the rate process and then, in Section 5, we introduced a generalization of lognormal processes to obtain the desired generalized fractional stochastic Gompertz curves.
A further extension to other growth curves with different rate functions could rely on nonlinear generalized fractional differential equations. While, on one hand, a theory of nonlinear generalized fractional differential equations is currently being developed (see, for instance [38]), to obtain stochastic growth curves of such type one could need some sort of nonlinear generalized fractional stochastic differential equations. In future works we will focus on generalizations of the Lévy-Liouville fractional Brownian motion (which is different from the well-known Mandelbrot-Van Ness fractional Brownian motion) to consider stochastic differential equations of the aforementioned type (see, for instance [23] for the fractional case). Let us also underline that the approach adopted in Section 4 can be easily extended not only to the multivariate setting, but also to stochastic processes defined on Banach spaces, the latter by using cylindrical noise processes, such as the cylindrical (possibly fractional) Brownian motion (see, for instance [51]). One can also define multivariate log Φ -normal distributions by a simple vectorization argument. However, once one moves from the one-dimensional to the n-dimensional setting, the discussion on the median cannot be reproduced as it is, but, instead, one should consider quantile contours (see, for instance [52]).
The models presented in this paper represent a natural step forward from [24,25]. Indeed, while in [24,25] the authors consider only Riesz kernels, here we propose a wide family of memory kernels (of which Riesz kernels are particular cases) that can be used. Moreover, analogously as done in [29] for the deterministic models, here one can consider two different memory kernels acting on the rate function and/or on the curve itself. Concerning the introduction of the noise, let us observe that stochastic models are shown to be useful, for instance, to study random fluctuations in mathematical oncology models (see, e.g., [8][9][10]). In [24] a fractional-order model has been used to describe tumour growth. Both the models presented in [29] and this paper provides some methods to introduce the noise in such fractional-order growth curves. There are three main differences between the present paper and [29]: • In [29] we considered only couples of Riesz kernels, while here we can consider any couple of suitable memory kernels; • Independently of the fact that memory effects are introduced in the rate function and/or in the curve, the generalized fractional stochastic Gompertz curve presented here is non-negative (as shown in Proposition 9), while this is not true in [29] in the case in which the curve function is defined itself via a fractional differential equation; • In [29], when the curve is defined via a fractional differential equation, the deterministic model is re-obtained by considering the mean of the stochastic curve. Here, in any case, the deterministic model is provided by the median of the stochastic curve, as in the classic case.
For these reasons, we think that the generalized fractional stochastic Gompertz curves defined in this paper are more realistic and more general models to describe Gompertztype growth phenomena with memory and noise, thus they represent a step forward with respect to [29]. However, restricting the view to a specific (parametrized) family of Bernstein functions permits a better calibration of the used model. Hence, before using the models we discussed here, it is advisable to consider an ansatz on the couple of Bernstein functions to consider.
In both deterministic and stochastic models, the tools to obtain such generalizations are provided by fractional and generalized fractional calculus, with particular attention to inversion formulae given in [36] and Grönwall-type inequalities (see [38,39]).
Finally, let us remark that the aim of the paper is to present a family of models that can be used for population growth phenomena with memory and/or noise. As evidenced by the cited literature, these kinds of phenomena naturally arise in physics, engineering, biology and social sciences and then more general models are useful tools for improving the knowledge about them.