A generalized Gompertz growth model with applications and related birth-death processes

In this paper, we propose a flexible growth model that constitutes a suitable generalization of the well-known Gompertz model. We perform an analysis of various features of interest, including a sensitivity analysis of the initial value and the three parameters of the model. We show that the considered model provides a good fit to some real datasets concerning the growth of the number of individuals infected during the COVID-19 outbreak, and software failure data. The goodness of fit is established on the ground of the ISRP metric and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d_2$$\end{document}d2-distance. We also analyze two time-inhomogeneous stochastic processes, namely a birth-death process and a birth process, whose means are equal to the proposed growth curve. In the first case we obtain the probability of ultimate extinction, being 0 an absorbing endpoint. We also deal with a threshold crossing problem both for the proposed growth curve and the corresponding birth process. A simulation procedure for the latter process is also exploited.


Introduction
Constructing growth curves describing dynamic evolutions is relevant to several applied fields. Indeed, growth dynamics are exhibited by a large number of real systems, such as the number of individual infected during an epidemic, the size of a living body during juvenility, other indicators associated with population growth, etc. The basic reference in this respect is the well-known Malthusian model, governed by the Revision submitted to Ricerche di Matematica on October 4, 2020.
B Antonio Di Crescenzo adicrescenzo@unisa.it Extended author information available on the last page of the article differential equation d N M (t) dt = r N M (t), t > 0, (1) where N M (t) represents the population size and r > 0 is the growth rate. This leads to exponential curve, adopted often in population ecology to describe growth in absence of constraints. Mathematical models associated to phenomena governed by growth tendency are based on suitable refinements of the exponential curve, such as the logistic and its generalizations treated in Tsoularis and Wallace [35]. Indeed, the traditional models are not always appropriate to describe certain growth phenomena since they involve few parameters. On the contrary, models involving more parameters are much more flexible and provide in general a better fit of the observed data. In this framework it is worth mentioning the contribution by Wu et al. [37], that adopted a model based on the generalized Richards curve for the COVID-19 outbreak. The latter curve is governed by the differential equation where N R (t) represents the cumulative number of cases at time t, r is the growth rate at the early stage, C is the carrying capacity (i.e., the final epidemic size), p ∈ [0, 1] and α are suitable shape parameters. Other forms of generalized growth models can be found in further investigations. For instance, Rincón et al. [29] analyzed a generalized Fujikawa's growth model described by the equation where r , α, γ , c, C and N m are constants, with r > 0, K > N m > 0, and C is the carrying capacity. Another case of interest has been treated recently by Chakraborty et al. [9] to propose a generalization of simple equations modeling the growth mechanism of biological processes, and finalized to generate more flexible shapes. Furthermore, a detailed treatment of extensions of the Gompertz-type equation in modern science has been presented in the book by Kyurkchiev and Iliev [24]. Along the lines of the above mentioned researches, the present paper is aimed to propose a suitable extension of the celebrated Gompertz model. This is a well-known growth model that is frequently adopted among the sigmoid models for fitting real data, and is governed by the following differential equation: with α, β > 0. Various re-parametrisations of this model have been considered by Tjørve and Tjørve [34]. The generalization proposed herewith is based on a suitable exponentiation of the term in parenthesis on the right-hand-side of (2). This leads to a more flexible growth model, which includes a variety of cases for the asymptotic analysis. Indeed, it may tend to infinity, or to a finite limit (the carrying capacity), or to zero. The resulting model provides also a generalization of a modified Korf model investigated recently in Di Crescenzo and Spina [12]. For a complete analysis of the model we study various features of interest, such as the correction factor, the relative growth rate, the inflection point, the maximum specific growth rate and the lag time. For the proposed model we also face a threshold crossing problem and perform a sensitivity analysis of the parameters. We purpose to analyze also a stochastic counterpart of the proposed model. This is motivated by the need of describing growth phenomena subject to random perturbations by means of mathematical models based on stochastic processes. Specifically, a desirable characteristics is that the mean evolution of the process coincides with the deterministic curve of the proposed model. In this framework, stochastic processes for the modeling of real growth phenomena have been largely considered in the literature. We recall the well-known approach based on diffusion processes for the stochastic model of tumor growth, such as that exploited in Albano and Giorno [1], Giorno et al. [18], Giorno and Nobile [16], Hanson and Tier [20], Spina et al. [31]. Other studies including Gompertz and logistic growth models based on stochastic diffusions can be found in Campillo et al. [8], Himadri Ghosh and Prajneshu [22], and Yoshioka et al. [38]. Recent advances involving fractional Gompertz growth models in biological contexts have been analyzed in Ascione and Pirozzi [4], Dewanji et al. [11], Frunzo et al. [15], and in Meoli et al. [25]. Nevertheless, our analysis will be restricted to the case of birth-death processes. The usefulness of this family of stochastic processes in the applications to sciences has been described recently in Crawford and Suchard [10]. We recall some previous investigations oriented to the analysis of birth-death processes for logistic and Gompertz stochastic growth, such as Di Crescenzo and Paraggio [13], Parthasarathy and Krishna Kumar [27], Swift [32] and Tan [33]. Various computational issues and a first-passage-time problem for time-inhomogeneous birth-death processes are investigated in Giorno and Nobile [17].
First we study an inhomogeneous birth-death process with linear time-dependent birth and death rates. The analysis is also developed towards the more interesting case of a time-inhomogeneous linear birth process. In both cases we specify the conditions that allow the mean of the process to be identical to the proposed generalized Gompertz growth curve.
In order to pinpoint the usefulness of the proposed model and the given results we focus on some applications to real data. Indeed, we show that the considered growth model is able to provide a good fit to various datasets in real cases of interest. As a first application we consider the data of the active and of the total number of contagions over a given period of the COVID-19 outbreak spread in Iran and Italy. Specifically, we analyze the proposed model, the Gompertz model and the logistic model as candidates to fit the considered data. We show that, among the cases under investigation, in many instances the proposed model provides the better fit under two comparison criteria. Indeed, in order to assess the goodness of the curve fitting, we adopt the ISRP metric (introduced recently in [5]) and the d 2 -distance. We point out that this investigation is finalized to perform a detailed investigation on an extended Gompertz growth model, rather than to solve the complex problem of finding a proper mathematical model to describe the evolution of the COVID-19 outbreak. The second application is devoted to software failure data from Tandem Computers. Also in this case, it is shown that the proposed model provides the better fit of the considered data under the ISRP metric and the d 2 -distance. This confirms that the model is effective for various instances of growth dynamics.
Let us now describe the plan of the paper. In Sect. 2 we introduce the model and illustrate its mean features, including the different shapes exhibited by the growth curve according to the parameters. Section 3 is devoted to a detailed analysis of the relevant features of the model. In Sect. 4 we then provide the announced applications of the proposed model to the number of outbreak contagions and to software failure data. In Sect. 5 we analyze a special inhomogeneous linear birth-death process, and provide a condition on the time dependent terms of the birth and death rates such that the mean of the process is equal to the proposed growth curve. A special case is treated in Sect. 6, where the case of a pure time-inhomogeneous birth process is considered. We also face the first-passage-time problem of this process through constant boundaries and time-varying boundaries. The latter case is developed by means of a simulation-based approach. To this aim a procedure able to simulate the birth times of the process is sketched. Some concluding remarks are finally provided in Sect. 7.

The proposed model
In this paper, we propose and study the growth model N G P D (t), which is solution of the following differential equation: with y > 0, A > 0, a > −1, a = 0 and b > 0. This is a suitable extension of the Gompertz model governed by Eq. (2). Moreover, such a model is included in the family of growth models N (t), t > 0, that are described by a different kind of differential equation, namely where N (t) denotes the size of a population at time t and ξ(t) > 0 is a time-dependent growth rate function. Suitable choices of the growth rate ξ(t) allow us to define a variety of models of interest. For instance, • the constant rate ξ(t) = r yields the classical Malthusian growth (cf. Eq. (1)): • the rate ξ(t) = ξ G (t) := αe −βt refers to the Gompertz growth model (cf. Gompertz [19]) • the rate ξ(t) = ξ K (t) := αt −(β+1) is concerning the Korf growth model (cf. Korf [23]) with α, β > 0, leads to a recently proposed growth model (cf. Di Crescenzo and Spina [12]) Let us now consider a suitable family of growth models, whose growth rate function is defined as follows: whereF(t) = P(X > t) is the survival function of a given nonnegative continuous random variable X , and A > 0 is a constant. From the differential equation (4) and Eq. (8), we have that the corresponding growth curve is given by Hereafter we show that the representation (9) allows to express some properties of N (t) in terms of the distribution of X .

Remark 1
The ultimate behaviour of the model (9) depends on the nature of the expectation of X , denoted as m := E[X ] = ∞ 0F (τ ) dτ . Indeed, one has where C is finite and denotes the carrying capacity of N (t).
Let us now recall the notion of increasing concave ordering (see Section 4.A of Shaked and Shanthikumar [30]). Given two random variables X i , i = 1, 2, having distribution functions F i (t) = 1 −F i (t), i = 1, 2, we say that X 1 is smaller than X 2 in the increasing concave order (denoted by Roughly speaking, this means that X 1 is both "smaller" and "more variable" than X 2 in some stochastic sense.

Remark 2
The increasing concave order allows to compare growth models of the form (9). For i = 1, 2, let It is not hard to see that if y 1 ≤ y 2 , A 1 ≤ A 2 , and X 1 ≤ icv X 2 , then N 1 (t) ≤ N 2 (t) for all t ≥ 0.
In order to propose a flexible growth model representing a generalization of various previous instances, we assume that X has a generalized Pareto distribution (GPD) with survival functionF where a > −1, a = 0, b > 0. We recall that some characterization results on the GPD can be found in Asadi et al. [3], in Section 4 of Hashemi et al. [21] and in Section 3.1 of Arriaza et al. [2]. In particular, this family includes three distributions, depending on the values of a: -the Pareto distribution, when a > 0; -the Power distribution, when −1 < a < 0 (in this case the distribution is bounded above); in particular, the distribution is Uniform if a = 1 2 ; -the Exponential distribution, when a → 0.
We point out that the curve (12) exhibits different behavior according to the value of a > −1, a = 0. In particular, by studying the derivative the following cases arise: (i) If a > 0, the curve (12) is well defined for all t ∈ (0, +∞); it is an increasing function for all t > 0; for t → +∞ it tends to the carrying capacity (ii) If −1 < a < 0 and 1 |a| is an odd integer, N G P D (t) defined in (12) is well defined for all t ∈ (0, +∞), with N G P D ( b |a| ) = C = y e Ab ; moreover, the curve N G P D (t) is increasing for all t > 0, and it goes to infinity for t → +∞. Fig. 1 The curve N G P D (t), given in (12), is plotted for y = 1, A = 1, b = 1 and a = −0.5 (solid), then it is decreasing for t ≥ b |a| and tends to zero for t → +∞. (iv) If −1 < a < 0 and 1 a is a non-integer real number, In Fig. 1, the curve (12) is plotted for different choices of a, with case (i) for a = 1, case (ii) for a = −0.2, case (iii) for a = −0.5, and case (iv) for a = −0.3. Note that when a = −0.3, for case (iv), the curve is defined only in 0 < t < 3.3 = b/|a| as specified above. Moreover, if a = −0.5, the random variable X considered in (8) is uniformly distributed; in this case the population size N G P D (t) tends to zero, i.e. to the extinction. (10) it is not hard to see that for the GPD distribution one has that t 0 F(τ ) dτ is increasing in a ∈ (−1, 0) ∪ (0, +∞) and is decreasing in b ∈ (0, +∞) for all t ≥ 0. Hence, denoting by N G P D,i (t) the growth model (12) characterized by parameters

Analysis of the growth model
In this section we discuss various features of the generalized Gompertz growth model proposed in (12).

The correction factor and the relative growth rate
In the context of growth analysis, interest is given to models governed by equations of the following form: The function f is a function of t only through the population size N (t), and is named a size covariate model. It allows to express the density dependent growth of the curve, and it is strictly related to the relative growth rate g (see Tsoularis and Wallace [35]) through the following relation: We remark that an extension of the relative growth rate, named modified relative growth rate, has been proposed and studied recently in Pal et al. [26] for the data analysis of growth models. Moreover, the function f is also called correction factor because it provides information about the deviation of a growth model from the classical exponential growth (for which f (z) = 1, z ≥ 0). From Eqs. (3) one has that for the GPD model the correction factor is given by As example, Fig. 2 shows some plots of the correction factor for model (12).

The inflection point
The monotonicity of the proposed model (12) has been studied in Sect. 2. Now we focus on the second derivative of N G P D (t) in order to compute the inflection point, when existing. From (13), the second derivative of N G P D (t), t > 0, results: For the study of the sign of (15), by considering the analysis shown in Sect. 2 we report two cases: -If Ab > a + 1, in cases (i), (iii) and (iv), the curve N G P D (t) is sigmoidal with the inflection point given by On the contrary, in case (ii) the curve N G P D (t) has two inflection points, i.e. t I and b |a| , with t I < b |a| . In particular, it has an upward concavity up to t = t I , a downward concavity between t I and b |a| , then an upward concavity. -If Ab < a + 1, in cases (i), (iii) and (iv), the curve N G P D (t) has a downward concavity for all t > 0, whereas in case (ii), N G P D (t) has a downward concavity up to the inflection point b |a| , then it has an upward concavity.
Furthermore, in all the interesting cases, the population at time t = t I is given by where C is given in (14).

The maximum specific growth rate and the lag time
The study of a growth curve that exhibits a sigmoidal behaviour in proximity of its inflection point (by approximating the curve in that point with a line) can be of interest in many phenomena that show lag, growth, and asymptotic phases. Due to the results of Sect. 3.2, we focus only on the instances of interest, that are cases (i), (iii) and (iv) analysed in Sect. 2, when Ab > a + 1. We consider the maximum specific growth rate, say μ, which is the coefficient of the tangent to the curve N G P D (t) given in (12) in the inflection point t I shown in (16): with μ expressed in (17). (For notation clarity, we denote by n the y-axis). Moreover, we introduce the lag time λ as the x-axis intercept of the tangent line. The lag time λ for the model (12) is the t-axis intercept of (18), that is Note that λ is positive when Ab > (a+1) 1 a +1 . Hence, the sigmoidal function N G P D (t) evolves with a growth rate that starts at zero and then accelerates to the maximal value μ, given in (17), in the time period resulting in the lag time λ shown in (19). Examples of tangent lines (dotted curves) and lag times (black point) are plotted in Fig. 3.

Threshold crossing
In several contexts of population evolution, it is relevant to know how long time the population spends below (or above) a certain preassigned value, which can represent a control threshold or a boundary in the evolution of the given phenomenon. Let us denote by S such a threshold for a growth model N (t). We are interested in the time instant in which N (t) reaches S, with S > N (0) = y. Taking into account the various behaviours that the curve N G P D (t) may have for different choices of a, we are interested in a generic threshold S > y if the curve grows to infinity, or on a percentage p of the carrying capacity C = ye Ab , or on a percentage of the size of the population at the inflection point t I given in (16). For a generic S, to obtain the crossing time instant θ S of N G P D (t) through S, we solve the equation N (θ S ) = S, recalling (12).
If S is a percentage 0 < p < 1 of the carrying capacity C = y e Ab , substituting S = pC in (20), the time of interest is Moreover, when S is a percentage 0 < p < 1 of N G P D (t I ), i.e. S = p y e Ab−a−1 , one has that the crossing time is: In Fig. 4, the quantities θ C ( p) and θ I ( p) are plotted as a function of p, with some choices of the parameters.

Sensitivity analysis
This section is devoted to analyze how the perturbation on each parameter involved in the model (12) influences the growth of N G P D (t). The parameters are the initial state y > 0, then A > 0, a > −1, a = 0, and b > 0. Recalling the study performed in Sect. 2 about the conditions for the existence, the following analysis can be performed on N G P D (t), that hereafter will be denoted by N ν G P D to emphasize the dependence on a generic parameter ν. Specifically, starting from (12) we expand N ν+ G P D in a Taylor series evaluated at ν, with > 0, for ν = y, A, a and b.
The latter term is positive for all t > 0.
To study the sign of (21), we note that due to the analysis performed in Sect. 2, it is not hard to see that The sign of the right-hand-side of (22) is equal to that of the function We have h(0) = 0. The first derivative h (t) is positive for all t > 0 in cases (i), (iii) and (iv) analyzed in Sect. 2; hence, in these cases one has h(t) > 0 for all |a| . Hence, the function h(t) is surely positive up to t = b |a| . Noting that lim t→+∞ h(t) = −∞, from the continuity of h we have that there exists at > b |a| such that h(t) < 0 for all t >t.
As example, in Fig. 5 we show the curve N G P D (t) and the effect of the perturbation = 0.1 on the parameter y (on the left) and on A (on the right). The same is show in    Fig. 6 for the parameter a, when a > 0 (on the left) and −1 < a < 0 (on the right), and finally in Fig. 7 for the parameter b in case (iii) (on the left) and in case (ii) (on the right).

Applications to real data
The usefulness of the results presented in the previous part of the paper emerges in various contexts in which the proposed generalization of the Gompertz model may be employed to obtain proper fit of real data. Indeed, in this section we show that the proposed GPD model (12) is able to fit well some datasets of interest. Two suitable criteria are then adopted in order to discuss the goodness of fit, by taking into account also other popular growth curves, i.e. the Gompertz model (5) and the logistic model The considered data analytic examples are related to epidemiology and reliability contexts. Indeed, the three considered models are used to fit the dataset related to (a) the active number and the total number of COVID-19 contagions from the 25th February to the 1st April 2020 in Iran and Italy (cf. [39] and [40], respectively); (b) a set of Release #1 failure data coming from software products at Tandem Computers (cf. Wood [36]), reported in Table 11. We use the nonlinear regression to fit the data, in particular the routine lsqcurvefit of Matlab ® , to solve nonlinear curve-fitting (data-fitting) problems in least-squares sense leading to the parameters estimation. Note that, we do not consider the Korf model (6) because for both datasets it does not perform a good fit, due to the fact that such a model has a slower growth when time increases. Moreover, we do not show the results related to the DS model (7) since it is a special case of the new model N G P D , and thus the fitting is better for the latter one.
Aiming to identify the best fitting model, we use two procedures to estimate parameters and check the goodness. In particular, we compute the sum of squared residuals (SSR) with the routine lsqcurvefit and then we estimate the parameters by means of ISRP metric. Further we calculate the d 2 -distance (in the euclidean norm) between ISRP and the constant parameter. We recall that the ISRP growth metric has been introduced recently in [5] aiming to determine the true growth curve that best fits the data (statistically). Indeed, such a metric provides an estimate of the rate parameter corresponding to the identified growth model in specific time intervals. The ISRP for the three considered models are respectively: with reference to the growth parameter α, r , A.
(a) We consider two data collections: (a1) consists of COVID-19 data n(i) collected every day from February 25th to April 1st, 2020, whereas (a2) is formed by the mobile mean of data collected in 3 days, i.e. data(i) = n(i)+n(i+1)+n(i+2) 35. The analysis of the mobile mean is suggested by the fact that the reliability of data on infected by COVID-19 is subject to high variability caused by external factors, such as fluctuations in the availability of test devices and in the test response times.
(a1) We show the interpolation of the dataset (a1), for three time intervals, where Figs. 8 and 9 refer to Italy and Figs. 10 and 11 refer to Iran. In all cases we have reasonably good fit. In Tables 1, 2 and 3, we report the SSR and the d 2 -distance between ISRP and the constant parameter, for the given datasets (a1). The notation e+05, for example, means ×10 5 . In Table 4, Total and Active Coronavirus Cases in Italy (March 29-April 1) are compared with the prediction done with the GPD model (12). The parameters used for the prediction at i-th day come from the estimation performed until up to the (i − 1)-th day, with the routine lsqcurvefit of Matlab ® . The same analysis is shown in Table 5 for Iran, whose Total and Active Coronavirus Cases (March 25-March 29) are compared with the prediction done with the GPD model (12). (a2) We show the interpolation of the dataset (a2), for three time intervals, for Italy in Figs. 12 and 13, and for Iran in Figs. 14 and 15. In all cases we have reasonably good fit. In Tables 6, 7 and 8, we report the SSR and the d 2 -distance between ISRP and the constant parameter, for the given dataset (a2). In Table 9, Total and Active Coronavirus Cases in Italy (March 28-March 31) are compared with the prediction done with the GPD model (12). The parameters used for the prediction at i-th day come from the estimation performed up to the (i −1)-th day. The same analysis is performed in Table 10, where Total and Active Coronavirus Cases in Iran (March 28-March 31) are compared with the prediction done with the GPD model (12).
We report in Tables 11 and 12 the estimates of the parameters for Figs. 8, 9, 10, 11, 12, 13, 14 and 15, and datasets (a1) and (a2), respectively, obtained by means of the routine lsqcurvefit of Matlab ® . By comparing the SS R and d 2 indexes of the two performed analysis, i.e. cases (a1) and (a2), we observe that (a2) provides better results. Indeed (a2) takes into account the running average of the data over 3 days. We find that according to the different intervals of time, and the adopted criteria (SS R and d 2 ), there is a variety in the detection of the best fitting model. Nevertheless, in the most of cases our proposed GPD model gives the best fitting. (b) We show in Table 13 the considered data, and in Fig. 16 we interpolate the data with the Gompertz, the Logistic and the GPD models. In Table 13 we report also the SSR and the d 2 -distance between ISRP and the constant parameters A, α, r , respectively. In these cases, we can observe that the GPD model fits better than the others.

Analysis of a special inhomogeneous linear birth-death process
Birth-death processes are largely used to model random evolution. See Callaert and Keilson [6,7] for the spectral structure of such processes. In this section, we consider a stochastic counterpart of the growth model introduced in (12) by studying an evolutionary model based on a birth-death process. We will consider a continuous-time Markov chain having an infinite state-space to mimics the growth curve (12) which can reach any high level when A is large. In particular, we consider a time-inhomogeneous linear birth-death process {X (t); t ≥ 0} having state space N 0 , and the absorbing endpoint 0. Denoting by the time-dependent transition rates of X (t), we assume that where λ(t) and μ(t) are positive functions, integrable on (0, t) for any finite t > 0. Note that the rates (24) are linear in i, hence the intensities of new births and deaths are proportional to the population size at the current time, and to the time-dependent functions λ(t) and μ(t), respectively. Clearly, λ(t) and μ(t) represent the individual birth rate and death rate at time t, respectively. We denote with P y,x (t) the transition probability of X (t), for all t ≥ 0, x ∈ N 0 and y ∈ N, i.e. the probability that the population, starting from y, reaches level x at time t. Note that, we consider X (0) = y ∈ N as for the growth model (12) the initial state is positive. Taking into account the results obtained in [12], the process X (t) with rates (24) has transition probabilities, for y ∈ N,  with m = min {y, x}, where ψ = ψ(t) and φ = φ(t) are given by: Moreover, the conditional mean E y (t) = E[X (t)|X (0) = y] and the conditional variance V y (t) = Var[X (t)|X (0) = y] of X (t) are respectively    Recalling the birth and death rates (24), the population mean satisfies the following differential equation: where ξ(t) is the net growth rate per capita of individuals, i.e.
Hence, by noting that equations (29) and (4) have the same form, we get that the mean of the process X (t) is equal to the growth curve proposed in (12), under the assumptions (11) and (30). Some properties of E y (t) and V y (t) are provided in Table  4 of [12]. Finally, we observe that P y,0 (t) is the probability that the population reaches the extinction prior to time t, being 0 an absorbing endpoint. In particular, from (25) the probability of ultimate extinction, conditional by the initial size y, is wherẽ (32) Obviously, ifμ −λ = ∞ orφ = ∞, then π y,0 = 1, i.e. the ultimate extinction is certain.

Analysis of a special case
In the following we assume that a > 0, that corresponds to the case (i) of the analysis performed in Sect. 2. This choice leads to a process X (t) having an increasing behavior and tending to a carrying capacity. As observed before, the conditional mean of X (t) verifies the same law of the growth model (12), under Eqs. (11) and (30). In this section, we focus our attention on this special case.

Proposition 1 The linear birth-death process X (t) with rates specified in (24) has conditional mean
if and only if Proof We substitute the first equation of (27) in (28) and we obtain

Hence, the expression (33) holds if and only if
that is equivalent to (34).
We now investigate on the process X (t) taking into account the relation (34). In this case, we have λ(t) > μ(t), and thus the net growth rate ξ(t) defined in (30) is decreasing and tends to 0 following a power law. Therefore, a population whose conditional mean has a S-shape is well described in this case. Moreover, for (27), the mean E y (t) is identical to the curve N (t) given in (12). Hence, due to (34), from (35) one obtainsλ −μ = Ab, withλ andμ defined in (32). From the results shown in Table 4 of [12], we have that E y (t) is strictly increasing, and it tends to the carrying capacity defined in (14), i.e. lim t→∞ E y (t) = ye Ab ≡ C. Moreover, we have withψ = e −Ab , so that V y (t) is strictly increasing in t.
Example 1 Under the validity of Eq. (34), we now consider various choices of the function μ(t) listed in Table 14. The transition probabilities (25) and (26) are plotted in Fig. 17 for y = 1, x = 0 and x = 1, and for some choices of the parameters. In Table 14 the last column is dedicated to the asymptotic absorption probability (31).
(a) We consider μ(t) = c, with c > 0, i.e. the individual death rate of the populations is constant. From (27) one has, for t > 0, where φ a (t) is defined in (37) and with k(t) = −Ab b at+b 1 a . The expression of the variance follows from (28).
where Q > 0 and c > |d| > 0. This case describes populations subject to individual sinusoidal death rate, due to i.e. periodic increase Table 14 Some choices of μ(t) and the corresponding extinction probabilities (31) Case , and the variance is obtained from (28).
In the first three cases of Table 14 one hasψ = +∞ becauseμ = +∞, therefore π y,0 = 1, i.e. the ultimate extinction is certain. On the contrary,μ < +∞ in case (d), and thusψ < +∞, so that π y,0 < 1. Moreover, we note that the variance diverges as t → ∞ for cases (a), (b) and (c), whereas it converges to a constant in case (d). In Fig. 18 we show some plots of V y (t), depending on μ(t).

Analysis of a special time-inhomogeneous linear birth process
In the previous section, we considered a birth-death process with conditional mean identical to the growth curve (12). Nevertheless the sample paths of that process may not reflect the behavior of the growth curve since they can be absorbed at zero (see the probabilities given in Eqs. (25) and (31)). In order to deal with a stochastic process more suitable to describe a growth phenomenon, we remove the possibility of downward jumps by assuming that μ(t) ≡ 0 in Eq. (24). This leads to a time-inhomogeneous linear birth process {X (t); t ≥ 0}, that possesses non-decreasing sample paths, and is characterized by time-dependent birth rates Here, S = {y, y + 1, . . .} is the state space and X (0) = y ∈ N is the initial state. Clearly, the function λ(t) is continuous, positive and integrable on (0, t), for any t > 0. As usual, we denote with P y, For y ∈ N and x ∈ S, one has (see, for example, [28]): where is the individual time-dependent cumulative birth intensity. By recalling Proposition 3 of Di Crescenzo and Spina [12] and the time-dependent growth rate (11), for a > 0, we have that the linear birth process with rates (38) has conditional mean if and only if This condition allows the conditional mean of the time-inhomogeneous birth process {X (t); t ≥ 0} to be equal to the growth curve given in Eq. (12). Moreover, from (40) and (42) we have and This property reflects that the births occur with decreasing time inputs, which is typical for environments with limited resources. Indeed, under condition a > 0 the process {X (t); t ≥ 0} reaches a mean saturation level. This is confirmed by the fact that the conditional mean (41) and the related variance are strictly increasing, with finite limits for Λ(t) given in (43). In Fig. 19 the mean (45) and the variance (46) are plotted for suitable choices of the parameters taken from the data considered in Sect. 4. We remark that the mean and variance exhibit a greater growth for the parameters considered in the time interval (a), that refers to a case in which the lockdown effects were not yet prevailing.
For this process, we also obtain the Fano factor (i.e. the index of dispersion, defined as the variance over the mean). Indeed, due to Eqs. (45) and (46), the following results hold: We thus obtain that D y (t) is monotonic increasing in t, with D y (0) = 0. Moreover, by analysing (47), we note that: is underdispersed, i.e. D y (t) < 1 for all t ≥ 0; in this case the occurrence of births is more regular than a Poisson process.
is underdispersed for t < t * , and overdispersed for t > t * , where therefore there is more irregularity in the distribution of births with respect to a Poisson process, for large times.
Now we consider the coefficient of variation of X (t), that from (45) and (46) results: where Λ(t) is given in (43); so σ y (t) is increasing in t, in A, a and b. Moreover, the following limiting behaviors are obtained: Therefore, when A → 0 or b → 0 there is a good correspondence between the deterministic law (12) and the stochastic process with birth rates (38). Clearly, if A → 0 or b → 0, then λ(t) → 0. This yields that for a low birth rate the two models exhibit a good agreement. It is worth mentioning that the knowledge of the transition probabilities (39) allows to perform a stochastic comparison for the birth processes considered in this section. To this aim we recall the following notion: Given two discrete random variables X i , i = 1, 2, we say that X 1 is smaller than X 2 in the likelihood ratio order (denoted by X 1 ≤ lr X 2 ) if P(X 2 = x)/P(X 1 = x) increases in x over the union of the supports of X 1 and X 2 (see Section 1.C of [30]). Then, in analogy with Remark 3 we are able to provide the following result. (10) it is not hard to see that for the GPD distribution one has that t 0 F(τ ) dτ is increasing in a ∈ (−1, 0) ∪ (0, +∞) and is decreasing in b ∈ (0, +∞) for all t ≥ 0. Hence, denoting by N G P D,i (t) the growth model (12) characterized by parameters y i , A i , a i , b i , i = 1, 2, due to Remark 2 we have that if y 1 ≤ y 2 , A 1 ≤ A 2 , a 1 ≥ a 2 and b 1 ≤ b 2 , then N G P D,1 (t) ≤ N G P D,2 (t) for all t ≥ 0.

First-passage-time problem
In Sect. 3.4, we considered the threshold crossing problem for the deterministic growth curve. Now, we analyze the similar first-passage-time problem of the considered stochastic process through certain thresholds. This is important to analyze relevant information on the reaching of critical levels in applied contexts. Let y ∈ S be the initial value and k ∈ N a threshold, with k > y. Aiming to analyze the first-passage-time random variable we denote by g y,k (t) = d P(T y,k ≤ t)/dt its probability density function (pdf). As in [12], we have where λ(t) and P y,k (t) are given in (42) and (39), respectively. For the behavior of g y,k (t) in proximity of t = 0 one has: lim t→0 g y,k (t) = A y , if k = y + 1 0, otherwise.
As example, in Fig. 20, the density (48) is plotted for some choices of the parameters. Moreover, Fig. 21 provides some instances of the mean of T y,k . For the first-passage-time problem of X (t) through a time-varying boundary we consider the continuous function t → β(t), where β(0) > y, and define T y,β = inf{t ≥ 0 : X (t) = β}, If β(t) is monotone nonincreasing, then similarly as in Proposition 3 of Di Crescenzo and Pellerey [14] we have with P y,n (t) given in (39). In more general cases the determination of the first-passagetime pdf g y,β (t) := d P(T y,β ≤ t)/dt can be handled by means of computational methods. In the next section we discuss a simulation procedure for X (t) which can be used to obtain estimates of first-passage-time pdf's.

Simulation
For the birth process {X (t); t ≥ 0}, characterized by time-dependent birth rates (38) and initial state X (0) = y ∈ N, we can construct a customary event-based simulation procedure. Denoting by T k the k-th increment (birth) of the process, with T 0 = 0, for all k ∈ N one has where Λ(t) is given in (43). This is not a bona fide distribution since, due to Eq. (44), Hence, given that T k = τ ≥ 0, the next birth occurs at a finite time t = T k+1 with probability 1 − , and thus the simulation of T k+1 can be performed through the steps indicated hereafter: 1. generate a random variable, say U , uniformly distributed in (0, 1); where, with A, a, b > 0, the inverse of Λ(t) is given by Figure 22 shows examples of simulated sample paths of X (t) obtained by means of the above sketched procedure. The latter may be used to develop a simulationbased approach to the first-passage-time problem of the birth process X (t) through a time-varying boundary β(t). Indeed, estimates of the pdf g y,β (t) can be constructed in terms of histograms obtained by means of extensive simulations based on the above procedure. As example, we consider the same case study treated in Sect. 4.1 of [14], for the boundaries (a) β(t) = log(t + 1) + 2 and (b) β(t) = 2 sin(π t/5) + 7. Estimates of the corresponding first-passage-time pdf's have been obtained through 10 5 simulated sample paths of X (t) performed by use of MATHEMATICA ® . See Fig. 23 for the corresponding histograms. As in similar cases exploited in the past, in case (a) the histogram exhibit changes of shapes when the boundary takes integer values, whereas in the case (b) the histogram reflects the periodicity of the periodic boundary. Moreover, we point out that the ratio of simulated sample paths of X (t) that reach the boundaries in the two cases is (a) 0.033 and (b) 0.202, respectively.

Concluding remarks
The choice of a suitable curve to describe a growth phenomenon is always a crucial task, since not all features of the relevant physical model and of the observed data can be captured by a given function. In many cases modelers are forced to perform the analysis on the ground of several curves aiming to compare the pertaining results and thus to detect the best choice on the basis of suitable statistical indexes. In some cases the presence of several parameters in the model allow to obtain a better fit of data, but an excess of parameters leads to a lack of correspondence between the considered model and the physical reality of the growth event. Hence, a proper compromise is required between the complexity of the model and its correspondence with the observed phenomenon. On the ground of these remarks, in this paper we proposed a growth model that provides a suitable generalization of the celebrated Gompertz model, but it is not excessively complex. Indeed, it involves three parameters, other than the initial (positive) value of the growth curve. After a thorough analysis of useful characteristics, we also focused on applications of the growth curve to real data concerning epidemiological and reliability contexts, where the proposed model is found very suitable to describe the observed dynamics under the ISRP metric and the d 2 -distance. We developed and analyzed two stochastic counterparts of the proposed model. They are based on an inhomogeneous linear birth-death process and a linear birth process. In both cases the correspondence between the growth curve and the mean of such stochastic processes is insured by a special choice of time-varying coefficients for the birth and death rates.