K R Philipsen

L E Christiansen

H Hasman

H Madsen

Modelling conjugation with stochastic differential equations

Keywords: SDE, horizontal gene transfer, resistance, CTSM, E. faecium, mathematical modelling, Hidden Markov Model, likelihood-ratio test, AIC

Conjugation is an important mechanism involved in the transfer of resistance between bacteria. In this article a stochastic differential equation based model consisting of a continuous time state equation and a discrete time measurement equation is introduced to model growth and conjugation of two Enterococcus faecium strains in a rich exhaustible media. The model

contains a new expression for a substrate dependent conjugation rate. A Maximum Likelihood based method is used to estimate the model parameters. Different models including different noise structure for the system and observations are compared using a likelihood-ratio test and Akaike's information criterion. Experiments indicating conjugation on the agar plates selecting for transconjugants motivates the introduction of an extended model, for which conjugation on the agar plate is described in the measurement equation. This model is compared to the model without plate conjugation. The modelling approach described in this article can be applied generally when modelling dynamical systems.

Introduction

Development and spread of antimicrobial resistance in bacterial populations is of increasing concern, as it can lead to major difficulties for the treatment of diseases. A first step in the direction of solving this problem is to gain a better understanding of the spread of resistance, and here conjugation plays a major role. Conjugation is one of several mechanisms of horizontal gene transfer by which plasmids coding for e.g. antimicrobial resistance can be transferred between bacteria. The use of mathematical modelling to describe the dynamics of plasmid spread and persistence was first introduced by [START_REF] Levin | The kinetics of conjugative plasmid transmission: fit of a simple mass action model[END_REF] and since then many studies have been made to improve the model framework and parameter estimation and to incorporate more accurate expressions for the plasmid dynamics [START_REF] Slater | Progress towards understanding the fate of plasmids in bacterial communities[END_REF]. Most studies have used ordinary differential equations, ODE [START_REF] Levin | The kinetics of conjugative plasmid transmission: fit of a simple mass action model[END_REF]Freter et al., 1983;[START_REF] Clewlow | Mathematical model of plasmid transfer between strains of streptomycetes in soil microcosms[END_REF][START_REF] Macdonald | Effects of energy availability on the conjugative-transfer kinetics of plasmid rp4[END_REF][START_REF] Willms | Static recipient cells as reser[END_REF], or ordinary difference equations [START_REF] Knudsen | Predictive model of conjugative plasmid transfer in the rhizosphere and phyllosphere[END_REF][START_REF] Sudarshana | Quantification and modeling of plasmid mobilization on seeds and roots[END_REF] as the modelling framework. This approach can be discussed since a part of the complexity of conjugation, e.g. dependence on the surrounding environment, is not included in these models. A way of overcoming this problem is to use a stochastic modelling approach where the randomness accounts for those processes not included in A c c e p t e d m a n u s c r i p t the model. Some efforts have been made to include randomness in plasmid models (Joyce et al., 2005;[START_REF] De Gelder | Stability of A c c e p t e d m a n u s c r i p t[END_REF][START_REF] Ponciano | The population biology of bacterial plasmids: A hidden markov model approach[END_REF].

Common for these studies is that they use discrete-time models, which is a good approximation at low bacteria densities, but not at higher densities where growth and conjugation is a continuous process. Joyce et al. (2005) gives a first and second order moment representation of their model which requires a derivation of the means and variances. Thus, it is not easy to transfer the method to other models. In the study by [START_REF] Ponciano | The population biology of bacterial plasmids: A hidden markov model approach[END_REF] one of the model parameters is implemented as a normal distributed random variable. This paper presents a generic modelling framework first described by Kristensen et al. (2004a), which is based on stochastic differential equations, SDE. The SDE is connected to data through a state-space formulation consisting of continuous-time state equations (the SDE) and discrete-time observation equations. The applied state-space approach is in fact a continuous time Hidden Markov Model. The state-space formulation opens up for strong statistical tools for estimating model parameters and for inference concerning the best model (Kristensen et al., 2004a). State-space modelling has also been used in previous plasmid studies [START_REF] De Gelder | Combining mathematical models and statistical methods to understand and predict the dynamics of antibiotic-sensitive mutants in a population of resistant bacteria during experimental evolution[END_REF], 2007;[START_REF] Ponciano | The population biology of bacterial plasmids: A hidden markov model approach[END_REF], but our approach differs by enabling continuous-time state equations in combination with discrete-time observations. The statespace model enables a simultaneous estimation of the growth and conjugation parameters, whereby the bacterial growth is accounted for when estimating the conjugation rate. This is an improvement from previous studies, where

A c c e p t e d m a n u s c r i p t

the donor and recipient concentrations were introduced as a mean value of two measurements taken over time [START_REF] Knudsen | Predictive model of conjugative plasmid transfer in the rhizosphere and phyllosphere[END_REF][START_REF] Sudarshana | Quantification and modeling of plasmid mobilization on seeds and roots[END_REF] or the experiment was constructed such that bacterial growth could be neglected [START_REF] Macdonald | Effects of energy availability on the conjugative-transfer kinetics of plasmid rp4[END_REF].

In this article the SDE based state-space modelling framework is applied to analyze data from an in vitro experiment for conjugation between two

Enterococcus faecium species growing in a batch exhaustable media. This experiment was made in order to study the transferability of vancomycin resistance in E. faecium. To describe the experimental system a conjugation model of bacteria growing in a broth exhaustible media is introduced, which is an expansion of previous models. Several authors [START_REF] Levin | The population biology of bacterial plasmids: A priori conditions for the existence of mobilizable nonconjugative factors[END_REF]Freter et al., 1983;[START_REF] Knudsen | Predictive model of conjugative plasmid transfer in the rhizosphere and phyllosphere[END_REF][START_REF] Clewlow | Mathematical model of plasmid transfer between strains of streptomycetes in soil microcosms[END_REF][START_REF] Top | Determination of the mechanism of retrotransfer by mechanistic mathematical modeling[END_REF][START_REF] Willms | Static recipient cells as reser[END_REF][START_REF] Sudarshana | Quantification and modeling of plasmid mobilization on seeds and roots[END_REF] have modelled conjugation events with the mass action model proposed by [START_REF] Levin | The kinetics of conjugative plasmid transmission: fit of a simple mass action model[END_REF]. The mass action model states that the appearance of transconjugants is proportional to the product of the donor and recipient concentrations. We

Experimental methods

A conjugation experiment was made in order to study the transferability of vancomycin resistance in Enterococcus faecium. As recipient the E.

faecium reference strain BM4105RF was used. This strain is resistant to rifampicin (MIC > 25 μg/ml) and fusidic acid (MIC > 25 μg/ml) due to chromosomally located mutations. As donor, the E. faecium A17sv1 (Hasman and Aarestrup, 2002) was used. This strain is resistant to erythromycin (MIC > 16 μg/ml) and vancomycin (MIC > 32 μg/ml) due to the presence of the erm(B) gene and Tn1546 transposon (carrying the vanA-gene cluster), respectively, located on a conjugative plasmid.

The conjugation experiment was performed in liquid Brain-Heart Infusion (BHI) media (Oxoid). Bacterial counting was performed on BHI agar plates supplemented with the following antibiotics when appropriate: rifampicin 25 μg/ml, fusidic acid (25 μg/ml), erythromycin (16 μg/ml) and vancomycin (32 μg/ml).

Conjugation experiment

Over night blood agar cultures of the two strains grown at 37 • C were inoculated in BHI media supplemented with the appropriate antibiotics as A c c e p t e d m a n u s c r i p t described above. From these, 100 μl of each culture was transferred to fresh tubes containing 10 ml preheated BHI media without antibiotics supplementation. When these cultures reached late exponential growth (OD 500 of 0.3-0.5), the number of cells in each culture was adjusted to the same amount and 1.5 ml of each was transferred to 100 ml preheated (37 • C) BHI media in an 250 ml Erlenmeyer flask. The flask was placed in a shaking incubator (125 rpm) at 37 • C. Four 1 ml samples were taken immediately after the cells were added (t=0) as well as 1. 50, 2.80, 3.00, 3.40, 3.80, 4.00, 4.40, 4.80, 5.00, 5.50, 6.00, 6.30 and 7.00 hours after the cells were mixed. One of the four 1 ml samples was used to measure the OD 500 of the culture. The remaining three samples were diluted ten-fold until 10 -7 . From here, 100 μl of each dilution (where appropriate) was plated onto sets of plates containing BHI agar supplemented with either of the following: 1) Recipient plates containing 25 μg/ml rifampicin + 25 μg/ml fusidic acid, 2) Donor plates containing 16 μg/ml erythromycin + 32 μg/ml vancomycin, and 3) transconjugant plates containing 25 μg/ml rifampicin + 25 μg/ml fusidic acid +16 μg/ml erythromycin + 32 μg/ml vancomycin. All plates were incubated 48 hours and the colony forming units, CFU, were counted. From these numbers the bacterial concentration in the flask was calculated.

Calculating bacterial concentration

The concentrations of the different bacteria populations in the flask can be determined from the plate count by a generalized linear model approach,

A c c e p t e d m a n u s c r i p t

assuming that the CFU count for each plate is Poisson distributed with an offset corresponding to the given dilution. If, for any given observation time t k , N jk is the count number j with dilution n j , then the expected CFU count E(N jk ) can be modelled with the generalized linear model

E(N jk ) = λ jk = n j exp(β k ) , N jk ∼ Poisson(λ jk ) ( 1 ) log(λ jk ) = log(n j ) + β k . (2)
Fitting the model to data gives an estimate for the coefficient β k and hereby an estimate of the bacterial concentration Y k = exp(β k ) in the flask at time, t k . The variance σ 2 Y,k of the estimated concentration at time t k can be determined from the variance of β as

σ 2 Y,k = Var[exp(β k )] = (exp(β k )Var[β]) 2 . ( 3
)
The concentrations Y k calculated in this way are the observations to be used for the modelling procedure described in the reminder of this article. The command glmfit in Matlab is used for fitting the generalized linear model (see Appendix A).

Model formulation

The modelling framework used in this study is a continuous-discrete time state-space model consisting of a continuous time state equation (the SDE) and a discrete time observation equation. The model has the general form A c c e p t e d m a n u s c r i p t (Kristensen et al., 2004b)

dX t = f (X t , u t , t, θ)dt + σ(u t , t, θ)dω t , (4)
Y k = h(X k , u k , t k , θ) + e k , (5)
where ( 4) is the SDE and ( 5) is the observation equation. X t is a ndimensional vector of state variables and Y k is a l-dimensional vector of observations. The observations are obtained at discrete times t k with the observation noise being e k ∈ N(0, Σ k ). θ is a vector of unknown parameters and u is a vector of input variables, i.e. variables which can be observed and have an influence on the system dynamics. The functions f and h can be linear as well as nonlinear functions. {ω t } is a standard Wiener process representing sources of noise in the system. The first part on the right-hand side of the SDE (4) is called the drift term and the second part is called the diffusion term. The stochastic state-space model ( 4)-( 5) has several advantages compared to deterministic models. For instance, the state-space approach separates the residual noise into system noise and observation noise, where the system noise is used to e.g. compensate for those biological processes into transconjugants with a rate dependent on the concentration of donor and recipient. In earlier studies [START_REF] Macdonald | Effects of energy availability on the conjugative-transfer kinetics of plasmid rp4[END_REF][START_REF] Knudsen | Predictive model of conjugative plasmid transfer in the rhizosphere and phyllosphere[END_REF][START_REF] Sudarshana | Quantification and modeling of plasmid mobilization on seeds and roots[END_REF] the substrate concentration and/or donor and recipient concentrations were kept constant in the experiments and were thus held constant in the model. In order to better approximate in vivo situations, these constraints were not applied to our E. faecium conjugation experiment, and thus all three populations as well as the substrate content must be included as state variables in the model. The drift term

A c c e p t e d m a n u s c r i p t

corresponding to the model shown in Figure 1 is

f (X t , θ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ μ D (S)D μ R (S)R -γ(S)DR μ T (S)T + γ(S)DR -η(μ D (S)D + μ R (S)R + μ T (S)T ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (6)
where μ i (S) is the growth rate for the bacteria population i (D, R or T ), and

γ(S) is the conjugation rate. X = [D, R, T, S] are the model state variables and θ = [v D , v R , v T , κ D , κ R , κ T , γ max , κ c , η]
are the parameters for the model.

The substrate, S, is simulated as a normalized variable with an initial value of 1. The amount of substrate in the solution decreases as it is utilized for bacterial growth. The parameter η is the amount of normalized substrate used for each cell division.

Bacterial growth continues until the substrate is exhausted, and the stationary state has been reached. The growth rate can be modelled with the well known Monod relation [START_REF] Monod | The growth of bacterial cultures[END_REF])

μ i (S) = v i S κ i + S , S ∈ [0, 1] , ( 7
)
where S is the substrate concentration, v i is the maximum growth rate and κ i is the substrate concentration when the growth rate is half of its maximum value. Both parameters are specific for a given bacteria population i.

The conjugation event depends on the probability of an encounter between

A c c e p t e d m a n u s c r i p t

donor and recipient. In the original mass action model the conjugation rate γ was introduced as a constant parameter. Several authors [START_REF] Knudsen | Predictive model of conjugative plasmid transfer in the rhizosphere and phyllosphere[END_REF][START_REF] Andrup | Kinetics of conjugative transfer: A study of the plasmid pxo16 from Bacillus thuringiensis subsp. israelensis[END_REF][START_REF] Andrup | A comparison of the kinetics of plasmid transfer in the conjugation systems encoded by the f plasmid from Escherichia coli and plasmid pcf10 from Enterococcus faecalis[END_REF][START_REF] Ponciano | The population biology of bacterial plasmids: A hidden markov model approach[END_REF] have discussed this assumption. [START_REF] Levin | The kinetics of conjugative plasmid transmission: fit of a simple mass action model[END_REF] showed that the mass action model presents a good estimation of the transconjugant population during exponential growth and under chemostatic conditions, but the model fails to describe the occurrences of transconjugants during the lag phase and at the onset of stationary phase. [START_REF] Macdonald | Effects of energy availability on the conjugative-transfer kinetics of plasmid rp4[END_REF] suggested that the conjugation rate depends on the substrate content and stated that conjugation can not occur without the presence of nutrition. However, to our knowledge this substrate dependence has not before been implemented in a mathematical model for conjugation. We suggest a nonlinear substrate dependent expression to model the conjugation rate

γ(S) = γ max S κ c + S , S ∈ [0, 1]. ( 8
)
This expression is similar to the Michaelis-Menten equation for enzyme kinetics and to the Monod relation. This expression is chosen as it forces the conjugation rate to reach a maximum value γ max when the substrate concentration is abundant, and it turns to zero as the substrate is depleted.

Depending on the value of κ c the conjugation rate will decrease concurrent with or after the decrease in the bacterial growth rate.

Some assumptions are applied to keep the conjugation model simple. (i)

A c c e p t e d m a n u s c r i p t

Transconjugants can function as donors transferring a plasmid to a recipient, but this is not described by the model. It is assumed reasonable to omit it for the E. faecium data set used in this study, as the concentration of transconjugants is very low compared to the donor concentration, and therefore does not contribute significantly to conjugation. (ii) The delay on 10-15 min that has been found [START_REF] Andrup | Kinetics of conjugative transfer: A study of the plasmid pxo16 from Bacillus thuringiensis subsp. israelensis[END_REF] between two conjugation events for the same donor is assumed to be insignificant. During the experiment there is always a large number of donors not involved in a conjugation event and thus ready to start the conjugation by encounter of a recipient. Therefore the delay is disregarded. (iii) It is assumed that the maximum growth rate of the transconjugants is either the same as the maximum growth rate for the recipient or smaller due to a fitness cost of the plasmid, i.e.

v T = v R (1 -α),
where α ∈ [0, 1[. An inference study will be made to test if α = 0.

The diffusion term

Depending on the system which is modelled it can be adequate to implement additive system noise, i.e. noise independent on the state variables, or multiplicative system noise, i.e. where the noise depends on the state variables. The choice of the system noise depends on assumption about the system modelled. [START_REF] Tier | Persistence in denisty dependent stochastic populations[END_REF] have described how different assumption of a biological process can lead to either demographic stochasticity (where the variance is proportional to the state variable) or environmental stochasticity (where the variance is proportional to the state variable A c c e p t e d m a n u s c r i p t squared). In our system the noise is implemented as multiplicative (environmental) noise, as the random fluctuations affect the whole population and not only the growth process. For instance the mass action model has shown to be good during exponential growth but not during lag-phase and stationary phase. Therefore multiplicative noise can be implemented to account for those processes not well described by the model. The method used for evaluating the likelihood function (Kristensen et al., 2004b) requires that the diffusion term is independent of the state variables. Therefore, instead of the state variables D, R and T , the input vectors u D , u R and u T , which contains observations for the three states, are inserted as a scaling for the standard deviation. The noise term is

σ m c = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ σ D u D 0 0 0 0 σ R u R 0 0 0 0 σ T u T 0 -ησ D u D -ησ R u R -ησ T u T σ S ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (9)
where an increment in donor, recipient or transconjugant concentration lower the substrate content in order to keep the mass balance, and σ S is introduced in order to ensure stability and is fixed to a small value.

An additive noise term σ a c is also introduced to compare with the multiplicative noise. This is done, as it is not sure that the data set is sufficient to make a good prediction of the system noise, and therefore a simplification of

A c c e p t e d m a n u s c r i p t

the noise structure can be an advantage. The additive noise term is

σ a c = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ σ D 0 0 0 0 σ R 0 0 0 0 σ T 0 -ησ D -ησ R -ησ T σ S ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 10
)

Observation equation

The observation equation relates observations to the state variables. In that way the state variables do not need to be measured directly and it is not necessary to have observations related to all state variables. In this study three of the state variables (donor, recipient and transconjugant) are observed. The concentration for the donor, recipient and transconjugant has been calculated from the CFU count as described in Section 2.2. We will refer to this as the observed concentration. The observation equation, when only conjugation in the broth media is considered, is

⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Y 1 Y 2 Y 3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ k = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ D R T ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ k + e k , e k ∈ N(0, Σ k ) . ( 11
)
The model with this observation equation is called the Broth model.

In addition to plasmid transfer in the broth, we have found that conjugation can occur on the transconjugant plates. This is a surprising result as antibi-

dR pk dt = -λR pk -γ p D pk R pk , (12)
dT pk dt = γ p D pk R pk . (13)
The initial values for the concentration of donor, D pk , and recipient, R pk , on the plate correspond to the concentrations in the flask at the time for the observation. Due to antibiotics on the plate, the donors and recipients die with a death rate λ and even though substrate is present in abundant amounts it is assumed that the donor and recipient bacteria can not grow.

The recipients which receives a plasmid on the agar plate (and thus become a transconjugants) become resistant to all the antibiotics present and can therefore grow and form colonies on the plate. The plate conjugation rate γ p is assumed to be independent off the substrate content as the substrate will not be depleted before the conjugation process stops.

The plate conjugation is modelled by finding the analytical solution for the

A c c e p t e d m a n u s c r i p t

ODE and inserting the solution in the observation equation for transconjugants Eq. ( 11). As the number of recipients receiving a plasmid on the plate is small, the last term in Eq. ( 13) can be disregarded. The equation to be solved is thus

dT pk dt = γ p D pk (t)R pk (t) , (15)
where

D pk (t) = D k0 exp(-λt) , (16)
R pk (t) = R k0 exp(-λt) . ( 17
)
The solution is

T pk = γ p 2λ D k0 R k0 (1 -exp(-2λt)) , (18)
which for t → ∞ is

T pk = γ p 2λ D k0 R k0 . (19)
The limit for t → ∞ can be considered since the CFU count is made after 12-24 hours, whereafter conjugation on the plate is very unlikely to occur due to death of the donor and recipients. Estimating both γ p and λ would be an over-parametrization of the model. Therefore γ p /λ is replaced by the

A c c e p t e d m a n u s c r i p t

parameter γ p . Inserting the result into Eq. ( 11) leads to the observation equation for the so-called Broth-plate model

⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Y 1 Y 2 Y 3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ k = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ D R T + γ p 2 DR ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ k + e k , e k ∈ N(0, Σ) . (20)

Observation noise

Three different covariance matrices for the observation equation are implemented and compared Additive:

Σ a k = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ s 2 1 0 0 0 s 2 2 0 0 0 s 2 3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (21)
Proportional:

Σ p k = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ s 2 4 σ 2 D,k 0 0 0 s 2 5 σ 2 R,k 0 0 0 s 2 6 σ 2 T,k ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (22)
Additive+Proportional:

Σ a,p k = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ s 2 1 + s 2 4 σ 2 D,k 0 0 0 s 2 2 + s 2 5 σ 2 R,k 0 0 0 s 2 3 + s 2 6 σ 2 T,k ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (23)

A c c e p t e d m a n u s c r i p t

Hence, the observations Y 1 , Y 2 and Y 3 are assumed to be uncorrelated in all cases. The variances σ 2 R,k , σ 2 R,k and σ 2 T,k are estimates from Eq. (3) at time t k .

The most simple noise form is the additive noise, but it is only reasonable, if the variance of the observations is independent of the observations. If this is not the case one alternative is to transform the data to stabilize the variance.

In this article instead a proportional noise term is suggested, which include the estimated variance of the observations. Additionally a covariance matrix is considered with both additive and proportional noise.

Statistical methods

The modelling procedure consist of several steps of parameter estimations and goodness of fit statistics. First the Broth model and Broth-plate models are reduced separately, i.e. the likelihood function is optimized for different implementation of the system noise, observation noise and drift term.

The best fit for the Broth and Broth-plate models are found applying goodness of fit statistics based on a ML approach. The inference study is made using a likelihood-ratio test and Akaike Information Criterion, AIC. After reducing the Broth and Broth-plate models they are compared, again using a likelihood-ratio test and AIC. Following the parameter estimation the models are simulated in Matlab, by implementing a numerical Euler method as described by Higham (2001).

A c c e p t e d m a n u s c r i p t

Parameter estimation

A ML estimation method is used to determine the SDE model parameters.

The parameters are found as those maximizing the likelihood function

L(θ; Y N ) = N k=1 p(Y k |Y k-1 , θ) p(Y 0 |θ) (24)
for the sequence of observations

Y N = [Y N , Y N -1 , . . ., Y 1 , Y 0 ]
. The conditional probability densities p are approximated by gaussian densities motivated by the fact that the SDE ( 4) is driven by a Wiener process having gaussian increments, i.e.

p(Y k |Y k-1 , θ) = exp -1 2 T k (Σ yy k|k-1 ) -1 k det Σ yy k|k-1 √ 2π l ( 25
)
where

Σ yy k|k-1 = Var[Y k |Y k-1 , θ] , ( 26
) k = Y k -Ŷ k|k-1 , and (27)
Ŷ k|k-1 = E[Y |Y k-1 , θ] . ( 28
)
The conditional mean Ŷ k|k-1 and covariance Σ yy k|k-1 in the likelihood function ( 24) and ( 25) can be estimated recursively by means of the Extended Kalman Filter (Kristensen et al., 2004b). In this study we use the software CTSM, from which the parameter estimation, correlation of the parameter estimates and the log-likelihood values are obtained. CTSM can be downloaded from the

A c c e p t e d m a n u s c r i p t

webpage: http://www2.imm.dtu.dk/∼ctsm/, from where a user manual is also available. CTSM is easy to use as it has a built in graphical user interface.

Recently an implementation of the method has also been made in R [START_REF] Klim | Population stochastic modelling (psm) -an r package for mixedeffects models based on stochastic differential equations, computer methods and programs in biomedicine[END_REF].

Goodness of fit statistics

Nested models are compared with a likelihood ratio test statistics given by

-2 log Λ = 2 ( (θ) -(θ 0 )) , ( 29
)
where the test statistics -2 log Λ is asymptotically χ 2 distributed with degrees of freedom equal to the difference in dimensions between the two models. (θ) = log(L(θ; Y N )) and (θ 0 ) = log(L(θ 0 ; Y N )) are the log-likelihood values of the model and the submodel, respectively. The inference study is also performed with AIC, which is given by

AIC = -2 (θ) + 2k , ( 30
)
where k is the number of parameters in the model and (θ) is the loglikelihood value of the model. When comparing models the preferred model is the one with the lowest AIC value.

A c c e p t e d m a n u s c r i p t

Results and discussion

The conjugation experiment continued for 7 h, during which samples were taken from the broth mixture approximately every 20 min and plated on selective agar plates. The CFUs were counted, and the bacterial concentrations in the broth were determined by a generalized linear model approach. The experimental results can be seen in Figure 2.

Inference study

Log-likelihood values for the Broth model and the Broth-plate model tested with respect to the state noise, the observation noise, and finally the drift term is summarized in Table 1. When estimating the model parameters only biological plausible parameter intervals are considered. The multiplicative system noise term (9) is tested against additive system noise (10), and it is found that there is no significant difference. We therefore continue the study for both the Broth and Broth-plate model with additive system noise. It would be interesting to investigate the system noise structure with a ML estimation based on the particle filter (Ionides et al., 2006), for which the multiplicative noise can be implemented directly dependent on the state variable. We will leave this for a future study. The additive system noise is well modelled with σ = σ D = σ R = σ T for both the Broth and Broth-plate models.

A c c e p t e d m a n u s c r i p t

Both the proportional observation noise and the additive observation noise can be tested against the proportional+additive observation noise using a likelihood-ratio test. The model with additive observation noise is seen to perform significantly worse than the proportional+additive noise model (pvalue = 0.0006 and p-value = 0.0007), whereas data is well modelled with

A c c e p t e d m a n u s c r i p t

proportional observation noise (p-value = 0.4467 and p-value = 0.4162). All the observation noise covariance matrices can be compared using AIC, from which it is also found that data is best modelled with proportional observation noise. This result is as expected, as the variance of the observations increases for higher CFU counts, which is best captured by proportional noise.

With this noise structure the two hypothesis α = 0 and κ T = κ R are tested.

We fail to reject both hypothesis as the p-values are high (between 0.7014 and 0.9969).

The best description of data for each of the models is thus a model with additive system noise where the standard deviation is the same for D, R and T .

The observation noise should be modelled as proportional noise and growth for the transconjugants and recipients can be modelled with the same set of parameters.

The reduced Broth and Broth-plate models can be compared with the likelihoodratio test, as they are nested models with the Broth-plate model containing one additional parameter (γ p ) compared to the Broth model.

Parameters

The result of the ML parameter estimation for the reduced Broth and Broth-plate models is shown in Table 2, and the results of simulations with these parameters are plotted in Figure 3 and4.

The conjugation rate is estimated to 2.008 • 10 -12 (CFU/μl) -1 h -1 for the Broth-plate model and 4.913 • 10 -12 (CFU/μl) -1 h -1 for the Broth model.

Thus, the choice of model influence the estimate of the conjugation rate.

The value of κ c (10 -16 -10 -15 ) is low compared to the κ value for the growth of the bacteria (10 -3 -10 -2 ). This means that the mass action model γDR 3. This means that the data is over-dispersed, and this over-dispersion is accounted for by σ 2 4 , σ 2 5 , and σ 2 6 . These parameters are in Poisson regression also refeered to as dispersion parameters. Several authors have found asymmetric likelihood profiles (Dennis et al., 2006;Ionides et al., 2006;[START_REF] Ponciano | The population biology of bacterial plasmids: A hidden markov model approach[END_REF][START_REF] King | Inapparent infections and cholera dynamics[END_REF] in state-space models for dynamical biological systems (e.g. bird growth, plasmid persistence and cholera pandemic).

In order to check the likelihood structure for our conjugation model, the profile-likelihood is calculated for v D and γ max . The profile likelihoods seen in Figure 5 asymmetric, which also in 3 out of 4 cases leads to asymmetric CI's. This is as expected due to the small number of observations. Furthermore it is seen that the profile likelihood CI is generally wider than the CI calculated by CTSM. However, this difference is not very large and we therefore believe that it has no or only a limited influence on the inference study performed. It should be noted that the observation noise seems to increase when the profile likelihood value is decreasing. This indicates that the observation noise (and not as expected the system noise) explains the difference between observations and the model. This is due to the few observations, which makes it difficult to adequately separate observation and system noise as discussed by Dennis et al. (2006).

In addition to the samples directly plated on the transconjugant plates, also 10 times concentrated samples were plated for the first six time points, where transconjugant concentration in the flask was low. However, an expected 100 times increase in conjugation on the plate compared to the non-concentrated sample was not observed. The reason for this is not clear. The sample is centrifuged in order to make a concentrated sample, this might change the ability for the donor and recipient to conjugate on the plate which could be one explanation. Further experiments should be performed to support the finding of conjugation under antibiotic pressure on the selective agar plate.

A c c e p t e d m a n u s c r i p t

  will introduce a new expression for the conjugation rate in an exhaustible media for which the proportionality constant of the mass action model is substrate dependent. An inference study is made to reduce the SDE model to its minimum form. A further extension of the model is made to treat adequately a methodological problem: the finding that conjugation continues on the agar plates selecting for transconjugants. In that case the observed transconjugants are a combination of transconjugants stemming from the the flask and conjugation occurring on the selective plates. The models with and without conjugation on the selective plates are compared in order to examine which model best describes data.

  not explicitly described by the model. The deterministic model often leads to autocorrelated residuals, which is not a problem with the SDE. Using the SDE based model also paves the way for strong statistical tools to estimate model parameters and make inferences.

Figure 1 :

 1 Figure 1: Flow diagram showing the structure of the conjugation models. The donor D can by encounter of a recipient R transfer a plasmid coding for resistance to the recipient. The recipient hereby is turned into a transconjugant T , which expresses the same antibiotic resistances as both the donor and the recipient.

  agar plate are traditionally considered to stop the conjugation process. However we have performed several experiments (data not shown), which have revealed that conjugation does also occur on the transconjugant plates. As the experiment is made to measure conjugation in the broth media we wish to separate conjugation on the plate from that in the flask. The following set of ODE's is suggested to model conjugation on the plate for each observation time k dD pk dt = -λD pk ,

Figure 2 :

 2 Figure 2: Data from the conjugation experiment with E. faecium growing in broth culture.

Figure 3 :

 3 Figure 3: The donor and recipient concentration over time as simulated by the Broth and Broth-plate models plotted together with data values.

Figure 4 :

 4 Figure 4: The transconjugant concentration plotted together with simulations of the Broth and Broth-plate models. For the Broth model the mean of the simulated transconjugant concentration and its confidence interval, as estimated by CTSM, is shown. For each model one example of the SDE simulation made in Matlab is plotted.

Figure 4

 4 Figure 4 we see that new transconjugants appear until approximately one hour after initiation of the stationary growth phase.When calculating the bacterial concentration from the CFU count it was assumed that the data was Poisson distributed, and thus that the variance equals the mean. If this was indeed the case each of the parameters σ 2 4 , σ 2 5 , and σ 2 6 should equal one. In this case σ 2 4 is around 14 and σ 2 5 and σ 2 6 are around

  are calculated by optimizing the likelihood function for fixed values of the parameter of interest. The 95% confidence interval (CI) plotted is the region of parameter values for which the profile log-likelihood value is larger than (θ) maxc/2, where (θ) max is the maximum log-likelihood value defined by Prob[χ 2 (1) < c] = 0.95. The profile likelihoods are indeed

Figure 5 :

 5 Figure 5: Profile likelihood (full line and circles) and approximate 95% confidence interval for v D and γ max . The cross marks the estimated parameter value. The dotted lines in the top plots give the standard deviation s 4 .

Table 1 :

 1 Log-likelihood values, p-values and AIC for different versions of the Broth model and Broth-plate model. The nested models listed below each other are tested with the likelihood-ratio test and the p-value for the comparison is listed in the line of the smallest of the two models tested. The proportional and additive observation noise structures are both tested against the proportional+additive noise term.

		Broth model	Broth-plate model
	( θ)	p-value AIC	( θ)	p-value AIC

Table 2 :

 2 The result of the ML parameter estimation for the Broth model and the Brothplate model. -16 (3.966 • 10 -16 ) 5.583•10 -15 (8.811•10 -15 ) γ max [(CFU/μl) -1 h -1 ] 4.913•10 -12 (0.636•10 -12 ) 2.008•10 -12 (0.699•10 -12 )

	Parameter	Broth model (SD)	Broth-plate	model
				(SD)
	v D [h -1 ]	1.678	(0.113)	1.787	(0.092)
	v R = v T [h -1 ]	1.216	(0.028)	1.261	(0.035)
	κ D	0.052	(0.047)	0.110	(0.047)
	κ R = κ T	0.0034	(0.0032)	0.0364	(0.0199)
	η	1.128•10 -6	(0.079 • 10 -6 ) 1.112•10 -6	(0.087•10 -6 )
	κ c 4.201•10 γ p [(CFU//μl) -1 ] 0 σ 4.010 • 10 -10 (2.828 • 10 -10 ) 1.128 • 10 -11 (2.488 • 10 -11 ) -4.096•10 -12 (1.630•10 -12 )
	σ 4	13.815	(5.349)	14.862	(5.575)
	σ 5	3.507	(1.212)	3.572	(1.304)
	σ 6	3.150	(1.422)	3.184	(1.280)

Acknowledgements

This study was supported by the Danish Research Council for Technology and Production Sciences through the grant 274-05-0117 titled "Evolution and adaptation of antimicrobial resistance in bacterial populations".


Patricia Román-Román

Desirée Romero

Francisco Torres-Ruiz

A diffusion process to model generalized von Bertalanffy growth patterns: Fitting to real data

Keywords: Diffusion process, stochastic growth, von Bertalanffy curve

The von Bertalanffy growth curve has been commonly used for modeling animal growth (particularly fish). Both deterministic and stochastic models exist in association with this curve, the latter allowing for the inclusion of fluctuations or disturbances that might exist in the system under consideration which are not always quantifiable or may even be unknown. This curve is mainly used for modeling the length variable whereas a generalized version, including a new parameter b ≥ 1, allows for modeling both length and weight for some animal species in both isometric (b = 3) and allometric (b = 3) situations.

In this paper a stochastic model related to the generalized von Bertalanffy growth curve is proposed. This model allows to investigate the time evolution of growth variables associated both with individual behaviors and mean population behavior. Also, with the purpose of fitting the above mentioned model to real data and so be able to forecast and analyze particular characteristics, we study the maximum likelihood estimation of the parameters of the model. In addition, and regarding the numerical problems posed by solving the likelihood equations, a strategy is developed for obtaining initial solutions for the usual numerical procedures. Such strategy is validated by means of simulated examples. Finally, an application to real data of mean weight of swordfish is presented.

A c c e p t e d m a n u s c r i p t 1 Introduction

Von Bertalanffy [START_REF] Bertalanffy | A quantitative theory of organic growth[END_REF] introduced an equation to study the growth of individuals belonging to several types of animal populations. As most growth models, it comes from an adaptation of the Verhulst logistic growth by assuming a maximal value of the growth variable (which might eventually be attained), and considering the growth rate as proportional to the difference between maximal and current value. It is currently the most common model used by fishery biologists to study growth in fish and its interpretations, such as fish population dynamics and the effects of fishery regulations on the catch.

The von Bertalanffy growth curve is

L(t) = L ∞ 1 -e -k(t-a) ; t ≥ a; k > 0, (1)
where L ∞ is the upper bound for the variable under study, that can only be reached after infinity time and k is the curvature parameter, or von Bertalanffy growth rate) that determines the speed with which the fish attains L ∞ . As regards the parameter a, sometimes named the initial condition parameter, it determines the time at which the fish has a size equal to zero and could be negative (in the biological literature is common to note this parameters as t 0 , but we have used a in order to avoid confusion with the notation usually employed for the initial instant in the context of stochastic processes). From a biological point of view, this question is meaningless because growth in the embryonic stage usually does not fit the von Bertalanffy growth pattern. Nevertheless, it is important to note that fish sufficiently aged to be exploited, as for consumption, show a trend modeled by the von Bertalanffy curve. The embryonic stage has then no interest in this context.

A general expression for the von Bertalanffy curve, also called "generalized von Bertalanffy growth curve" (see García-Rodríguez et al. [START_REF] García-Rodríguez | On the biology growth of the anglerfish lophius Budegasa Spinola, 1807 in the Spanish Mediterranean: a preliminary approach[END_REF] and references therein), is

B(t) = B ∞ 1 -e -k(t-a) b ; t ≥ a; k > 0; b ≥ 1, (2)
where the parameter b can be known or unknown. For example, the value b = 1 (curve [START_REF] Albano | A stochastic model in tumor growth[END_REF]) is used when the variable under study is the length. However when focusing on the weight, and taking into account the existing relation between the weight and the length, the value b = 3 is associated with the isometric growth, whereas the case b = 3 is related to the allometric growth.

Curves [START_REF] Albano | A stochastic model in tumor growth[END_REF] and [START_REF] Capocelli | Growth with regulation in random environment[END_REF] provide appropriate deterministic models to describe the growth (in length or weight) of fish and others animals. The study of procedures in order to determine the parameters of such curves for fitting real data

A c c e p t e d m a n u s c r i p t

has been widely considered. For example, Rafail [START_REF] Rafail | A simple and precise method for fitting a Von Bertalannfy growth curve[END_REF] developed a procedure based on a straight-line relationship between the natural logarithms of growth increments per unit of age against age as the independent variable.

Nevertheless, such models do not include variability among individuals of the same age, or environmental variability (fluctuations or disturbances that might exist in the system under consideration).

In order to take these variations into account, several stochastic models have been considered in the literature. The first type of model, according to the classification provided by Russo et al. [START_REF] Russo | Lévy processes and stochastic von Bertalanffy models of growth, with appliccation to fish population analysis[END_REF], are those described by a stochastic differential equation, obtained by including a noise term in the ordinary differential equation of the respective deterministic model.

In this context, and as concerns the Bertalanffy growth model, Lv and Picthford [START_REF] Lv | Stochastic Von Bertalanffy models, with applications to fish recruitment[END_REF] have recently considered three stochastic models, associated with the curve [START_REF] Albano | A stochastic model in tumor growth[END_REF]. Said models are built from stochastic differential equations with identical von Bertalanffy deterministic parts and different stochastic terms, to show, among other points, that stochasticity can have positive impact on fish recruitment. Gudmundsson's [START_REF] Gudmundsson | Stochastic growth[END_REF] criticism of such models is that the solution for these differential equations may not be strictly increasing. From this point of view, the aforementioned models may not be appropriate for the study of some growth variables, like length, that are not subject to decreasing (Weatherly and Gill [START_REF] Weatherly | The Biology of Fish Growth[END_REF]).

Such defect, real from a theoretical point of view, is practically diluted in practice. In the practical application of this type of models, the noise term has a moderate magnitude (see Gutiérrez et al [START_REF] Gutiérrez | A new gompertz-type diffusion process with application to random growth[END_REF] in the context of Gompertz-type growth). This feature makes the sample paths of the resulting processes show von Bertalanffy-like behavior patterns, except for the presence of small random fluctuations. In addition, if the resulting stochastic model verifies that its mean function is a von Bertalanffy growth curve, the use of the model for fitting and predicting is fully justified. On the other hand, such defect disappears if the stochastic model is associated with the von Bertalanffy generalized curve. This allows for the study of weight-to-age (generalizing the case for length-to-age) since the evolution of the animal's weight in time is not necessarily a strictly increasing function.

The second type of stochastic models are those that assume that the parameters of the von Bertalanffy curve are different for each member of the population and are thus considered random variables with a certain probability distribution.

Among the authors that deal with this second type of models, we must emphasize the one by Cheng and Kuk [START_REF] Cheng | Determination of the unknown age at first capture of western rock lobsters (panulirus cygnus) by random effects model[END_REF], that considered the parameters of the model as random effects following a trivariate normal distribution. More recently,

A c c e p t e d m a n u s c r i p t

Tóvar-Ávila et al. [START_REF] Tovar-Ávila | Use of stochastic models to estimate the growth of the Port-Jackson shark, Heterodontus portusjacksoni, off eastern Vistoria, Australia[END_REF] considered a reparametrization of the von Bertalanffy growth rate using three different probability distribution functions: Weibull, gamma and log-normal. Such models show an increasing behavior and consider individual variations, but are limited by the fact that they do not account for environmental variations, caused by multiple factors which are not always quantifiable or may even be unknown.

Russo et al. [START_REF] Russo | Lévy processes and stochastic von Bertalanffy models of growth, with appliccation to fish population analysis[END_REF], in order to overcome the flaws of the two types of models already mentioned, introduced a type of stochastic models related to the von Bertalanffy curve as a solution for stochastic equations that include a subordinator (a class of strictly increasing stochastic processes which make the solution process of the stochastic equation also increasing). Such models account for both individual and environmental sources of randomness.

In practical situations, the use of a stochastic model for fitting, making forecasts or drawing conclusions about the particular growth of the individuals of a population or about the mean population growth, requires the estimation of the unknown parameters of the model.

As a first approach along this line, some authors, such as Kimura [START_REF] Kimura | Likelihood methods for the Von Bertalanffy growth curve[END_REF], have studied likelihood methods under the assumption of independent and normally distributed errors using classic nonlinear least square methods. Later, Wang [START_REF] Wang | Estimating equations for parameters in stochastic growth model from tag-recapture data[END_REF] introduced unbiased estimating functions for a class of growth models that incorporate stochastic components and explanatory variables. More recently, Hart and Chute [START_REF] Hart | Estimating von bertalanffy growth parameters from growth increment data using a linear mixed-effects model, with an application to the sea scallop Placopecten magellanicus[END_REF] introduced a novel linear mixed-effects method for estimating von Bertalanffy growth parameters from growth increment data that lack explicit age information. In Cheng and Kuk [START_REF] Cheng | Determination of the unknown age at first capture of western rock lobsters (panulirus cygnus) by random effects model[END_REF] the problem of estimating the parameters of their already mentioned model is considered (see also Laslett et al. [START_REF] Laslett | Biological and Statistical Issues in Fitting Growth Curves to Capture-Recapture Data[END_REF] for a discussion of the results).

Russo et al [START_REF] Russo | Lévy processes and stochastic von Bertalanffy models of growth, with appliccation to fish population analysis[END_REF] also dealt with the problem of parameter estimation based on the probability distributions of the growth variable in each of the observed time instants.

Nevertheless, in the case of models described by stochastic processes that are a solution of stochastic differential equations, an efficient estimation of parameters must be based on data that provide information of the evolution of variables along time (sample paths). Thus such estimation must be based on data which can be of two types: In the first case, and using several sample paths, it is possible to estimate a model related to the particular growth of individuals, whereas in the second case, and from a single sample path, the estimated model is associated with the mean population growth.

Thus, the aim of this article is twofold. Firstly, the introduction of a stochastic model (specifically a diffusion process) to model behavior patterns associated with the generalized von Bertalanffy curve (2) that allows the consideration of both length and weight for some animal species. To this end, the model is obtained by applying the methodology developed by Gutiérrez et al. [START_REF] Gutiérrez | A new gompertz-type diffusion process with application to random growth[END_REF] in the context of gompertzian growth: from the deterministic equation, whose solution is the curve of interest, a stochastic component is introduced as well as the condition that the mean function of the resulting diffusion process be a curve of the type (2) (that fits well the sample data and can be used for forecasting purposes) is imposed. The building of the model as well as the study of some of its characteristics (probability distribution of the process, mean, mode and quantile functions) is presented in sections 2 and 3. For point predictions the mean and the mode functions can be used, whereas for interval predictions, the quantile functions provide intervals containing the growth variable of the process, for each time, with a specific probability. Furthermore, the consideration of a diffusion process in this context will allow the study of time variables, such as first-passage-times, associated with von Bertalanffy type growth models.

Secondly, since the purpose of model is to fit and forecast real data, the next objective is the estimation of the model. In section 4 an inferential study of the parameters of the process is carried out (particularly their maximum likelihood estimation) on the basis of discrete sampling (in both cases, by considering one or several sample paths). In addition, and regarding the numerical problems posed by solving the likelihood equations, in subsection 4.1 a strategy is developed for obtaining initial solutions in order to apply the usual numerical procedures. Such strategy is validated by means of simulated examples in section 5. Finally, in section 6 the possibilities of the new stochastic model, relative to fit, forecasting and first-passage-times, are illustrated by means of an application to real data of mean weight of swordfish in Southeastern Pacific.

A c c e p t e d m a n u s c r i p t 2 The generalized von Bertalanffy diffusion process

Before introducing the diffusion process, we will obtain an equivalent expression for the generalized von Bertalanffy curve (2), more suitable for our purpose. We will consider that the observations of the variable under study are made from a time instant t 0 ≥ 0, with x 0 > 0 the observed value at t 0 . With this realistic hypothesis in mind, and taking into account (2), we deduce that a < t 0 , since this function is increasing and verifies B(a) = 0. On the other hand, as we suppose B(t 0 ) = x 0 , and denoting c = e ka , we conclude that

B ∞ = x 0 (1 -ce -kt 0 ) b , (3)
from where

B(t) = x 0 1 -ce -kt 1 -ce -kt 0 b , t ≥ t 0 > ln c k ; k > 0; b ≥ 1. ( 4
)
Now we introduce the new diffusion process related to the curve (4). To this end, we follow the methodology considered in Gutiérrez et al. [START_REF] Gutiérrez | A new gompertz-type diffusion process with application to random growth[END_REF] in which several methods are presented in order to introduce a diffusion process in the context of gompertzian growth. Specifically, we look for a process in which the solution of the Fokker-Planck equation, without noise, is such a curve. In addition, the resulting process must verify that the mean function conditioned on the initial value, E[X(t)|X(t 0 ) = x 0 ], coincides with (4), which is specially useful for the purposes of forecasting.

Following the same scheme presented there, let us consider the first order equation verified by the transition probability density of the process, f ≡ f (x, t|x 0 , t 0 ),

∂f ∂t = - bck e kt -c ∂ ∂x [xf ] , x > 0, t ≥ t 0 with the initial condition lim t→t 0 f (x, t|x 0 , t 0 ) = δ(x -x 0 ). Its solution is f (x, t|x 0 , t 0 ) = δ ⎛ ⎝ x -x 0 1 -ce -kt 1 -ce -kt 0 b ⎞ ⎠ ,
which implies that the population under consideration increases according to [START_REF] Chong | Edad y crecimiento del pez espada (Xiphias gladius Linnaeus, 1758) en el Pacífico suroriental (Diciembre 1994-Septiembre 1996)[END_REF]. Now, we consider the Fokker-Planck equation of the homogeneous

A c c e p t e d m a n u s c r i p t lognormal diffusion process with infinitesimal moments

A 1 (x) = mx and A 2 (x) = σ 2 x 2 (m, σ > 0), ∂f ∂t = - ∂ ∂x [mxf ] + σ 2 2 ∂ 2 ∂x 2 [x 2 f ], x > 0, t ≥ t 0 ,
and modify the infinitesimal mean by multiplying it by the term bck/ e ktc . Thus, we obtain

∂f ∂t = - bck e kt -c ∂ ∂x [xf ] + σ 2 2 ∂ 2 ∂x 2 [x 2 f ], x > 0, t ≥ t 0 , (5)
that is, the forward equation of a new diffusion process with infinitesimal moments

A 1 (x, t) = bck e kt -c x (6) A 2 (x, t) = σ 2 x 2 .
It is obvious that the solution (5), when σ 2 vanishes, is the curve (4).

Alternatively, the process can be obtained from the Langevin equation

dX(t) dt = bck e kt -c X(t) + X(t)σW (t), (7)
where W (t) denotes the standard Wiener process The derivation of ( 7) can be achieved from the deterministic growth equation

dx(t) dt = bck e kt -c x(t), x(t 0 ) = x 0 , (8)
which can be seen as a generalization of the malthusian growth model with a deterministic fertility depending on the time, r(t) = bck/ e ktc , and replacing this fertility with r(t) + σW (t). Moreover, equation ( 7) has an important meaning in terms of how the population increases. In fact, this equation can be expressed in the form

dX(t) dt = mX(t)[1 -ϕ(X(t), t)],
where ϕ(x, t) = 1bck/ e ktc -(1/m) σW (t) is the well-known regulation function (see Capocelli and Ricciardi [START_REF] Capocelli | Growth with regulation in random environment[END_REF]) that, over time, introduces several changes into m, the growth rate of the Malthusian model.

A c c e p t e d m a n u s c r i p t

By rewriting [START_REF] Gutiérrez | First-passagetime densities for time-non-homogeneous diffusion processes[END_REF] in the usual form for stochastic differential equations, that is

dX(t) = bck e kt -c X(t)dt + σX(t)dW (t), (9)
its solution is a non-homogeneous diffusion process {X(t); t t 0 } taking values on R + and with infinitesimal moments

A 1 (x, t) = h(t)x (10) A 2 (x, t) = σ 2 x 2 ,
where h(t) = bck/ e ktc or h(t) = bck/ e ktc + σ 2 /2 according to whether the Itô or the Stratonovich integral is used to solve it, respectively. In addition, it is not difficult to show that the mean function, conditioned on the initial value x 0 , is

E [X(t)|X(t 0 ) = x 0 ] = x 0 exp ⎛ ⎝ t t 0 h(s)ds ⎞ ⎠ = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ x 0 1 -ce -kt 1 -ce -kt 0 b
for the Itô solution

x 0 1 -ce -kt 1 -ce -kt 0 b exp σ 2 2 (t -t 0 )
for the Stratonovich solution.

(

From [START_REF] Kimura | Likelihood methods for the Von Bertalanffy growth curve[END_REF], and taking into account our objectives, the usefulness of the model to fit and forecasting purposes, the Itô solution must be chosen (note that in the case of the Stratonovich solution, the mean is not even a bounded curve), so we introduce a new von Bertalanffy-type diffusion process associated with the curve (4) as the diffusion process {X(t); t t 0 } defined on R + and with infinitesimal moments given by [START_REF] Gudmundsson | Stochastic growth[END_REF].

Finally, another way to develop this process, in the line shown by Albano and Giorno [START_REF] Albano | A stochastic model in tumor growth[END_REF] in the Gompertzian case, is based on the discretization of ( 8) and then its randomization (see Appendix A for more details).

In short, the diffusion process presented is a stochastic model associated with the generalized von Bertalanffy curve (2), or its rewriting as (4), which accounts for both individual and environmental variability. Its main advantages are:

• Unlike other existing stochastic models, related to the von Bertalanffy curve (1) and useful to study the length of some animal species, the model under discussion allows to study both length and weight.

A c c e p t e d m a n u s c r i p t

• Its mean function is a generalized von Bertalanffy curve, which justifies it use for the fitting and prediction of real data with such a behavior pattern.

In addition, one of the advantages of using this sort of model for the study of dynamic phenomena is that it allows to consider some questions about their evolution through time variables. For instance, the study of time variables as the time an animal takes to reach the minimum size at which it can be sold for consumption, or the time a population takes to reach a certain size, can be of great interest. Problems like these may be solved by obtaining the density function of the first-passage-time through a constant boundary. The application developed in section 6 deals with a problem of this type, illustrating the model's usage for this end.

Probability distribution and some characteristics of the process

The probability distribution of the process, determined by the finite-dimensional distributions, can be obtained from the theory of stochastic differential equations, as solution of ( 9), or employing the theory of partial differential equations from the forward equation, [START_REF] García-Rodríguez | On the biology growth of the anglerfish lophius Budegasa Spinola, 1807 in the Spanish Mediterranean: a preliminary approach[END_REF], and the Kolmogorov or backward equation

∂f ∂t 0 + bck e kt 0 -c x 0 ∂f ∂x 0 + σ 2 2 x 2 0 ∂ 2 f ∂x 2 0 = 0, x 0 > 0, t ≥ t 0 .
In [START_REF] Gutiérrez | Approximate and generalized confidence bands for some parametric functions of the lognormal diffusion process with exogenous factors[END_REF] both approaches have been developed in the context of the lognormal diffusion process with exogenous factors, of which the von Bertalanffy process is a particular case. In this case, the transition probability density function of the process is

f (x, t|y, s) = 1 x 2πσ 2 (t -s) × exp ⎛ ⎜ ⎝- 1 2 ln x y -b ln 1-ce -kt 1-ce -ks + σ 2 2 (t -s) 2 (t -s)σ 2 ⎞ ⎟ ⎠ , t > s,( 12
)
which corresponds to a lognormal distribution, that is,

[X(t)|X(s) = y] ∼ Λ ln y + b ln 1 -ce -kt 1 -ce -ks - σ 2 2 (t -s); (t -s)σ 2 .( 13
)
Because of the Markovian property of the process, from ( 13) and the initial distribution, we can calculate the finite-dimensional distributions. In this

A c c e p t e d m a n u s c r i p t

case, we consider two initial distributions: a degenerate distribution, that is, P [X(t 0 ) = x 0 ] = 1, and a lognormal distribution, X(t 0 ) ∼ Λ(μ 0 , σ 2 0 ), these choices ensuring that the finite-dimensional distributions are lognormal. We must remark that the former choice can be seen as a particular case of the second considering σ 0 = 0 and μ 0 = ln(x 0 ). Moreover, the degenerate initial distribution is the real situation when only a sample path is available, whereas the lognormal case requires several trajectories. In any case, the random vector (X(t 1 ), . . . , X(t n )) follows an n-dimensional lognormal distribution Λ n (μ, Σ), where

μ i = μ 0 + b ln 1 -ce -kt i 1 -ce -kt 0 - σ 2 2 (t i -t 0 ), i = 1, . . . , n,
and

Σ ij = σ 2 0 + σ 2 (Min{t i , t j } -t 0 ), i, j = 1, . . . , n.
We now describe the main characteristics of the process, focussing particularly on the three most commonly employed in practice, especially for forecasting purposes. These characteristics are the mean function (which by the structure of the model is a type (2) curve, and thus particularly appropriate for fitting and predicting), the mode function (that provides, for each time, the most probable value of the growth variable) and the quantile functions (which allow to make predictions through intervals that contain the growth variable for each time, with a specific probability). Its expressions can be formulated jointly for the two initial distributions under consideration:

• Mean function m(t) = E[X(t)] = E[X(t 0 )] 1-ce -kt 1-ce -kt 0 b , t ≥ t 0 . • Mode function Mo(t) = Mode[X(t)] = Mode[X(t 0 )] 1-ce -kt 1-ce -kt 0 b exp -3 2 σ 2 (t -t 0 ) , t ≥ t 0 . • Quantile function C α (t) = α -quantile[X(t)] = α -quantile[X(t 0 )] 1 -ce -kt 1 -ce -kt 0 b × exp - σ 2 2 (t -t 0 ) + z 1-α σ 2 (t -t 0 ) + V ar[ln(X(t 0 ))] -V ar[ln(X(t 0 ))] , t ≥ t 0 ,
where z 1-α is the α-th quantile of a standard normal distribution.

A c c e p t e d m a n u s c r i p t

Furthermore, and from ( 13), the conditioned versions of these functions can be also obtained (see [START_REF] Gutiérrez | Approximate and generalized confidence bands for some parametric functions of the lognormal diffusion process with exogenous factors[END_REF] in the general case).

• Conditional Mean function m(t|s) = E[X(t)|X(s) = x s ] = x s 1-ce -kt 1-ce -ks b , t > s. • Conditional Mode function Mo(t|s) = Mode[X(t)|X(s) = x s ] = x s 1-ce -kt 1-ce -ks b exp -3 2 σ 2 (t -s) , t > s. • Conditional Quantile function C α (t|s) = α -quantile[X(t)|X(s) = x s ] = x s 1 -ce -kt 1 -ce -ks b × exp - σ 2 2 (t -s) + z 1-α σ 2 (t -s) , t > s.

Inference on the model

As we have noted above, the mean function of the process is a von Bertalanffy curve of the type (4). Therefore the mean function, as well as its conditional version, E[X(t)|X(s) = x s ], t > s, can be useful for making predictions with this model. For this reason, let us examine in this section its maximum likelihood (ML) estimation. First, we obtain the ML estimators of the parameters of the model and then that of the mean function, as well as the other parametric functions abovementioned.

Let us consider a discrete sampling of the process, based on d sample paths, for times t ij , (i = 1, . . . , d, j = 1, . . . , n i ) with t i1 = t 1 , i = 1, . . . , d. That is, we observe variables X(t ij ), the values of which, {x ij } i=1,...,d;j=1,...,n i , make up the sample for the inferential study.

The likelihood function depends on the choice of the initial distribution. When

P [X(t 1 ) = x 1 ] = 1, this function is L x ij (b, c, k, σ 2 ) = d i=1 n i j=2 f (x ij , t ij |x i,j-1 , t i,j-1 ),
where b, c, k and σ 2 are the parameters to be estimated. If

X(t 1 ) ∼ Λ(μ 1 , σ 2 1 ) the likelihood is L x ij (μ 1 , σ 2 1 , b, c, k, σ 2 ) = d i=1 f X(t 1 ) (x i1 ) n i j=2 f (x ij , t ij |x i,j-1 , t i,j-1 ).
In the second case there are two additional parameters that must be included in the estimation procedure. Nevertheless, the estimations of μ 1 and σ 2 1 depend only on the initial values and do not influence the estimation of the other parameters. Hence, the ML estimators of b, c, k and σ 2 are the same in both cases.

A c c e p t e d m a n u s c r i p t

Henceforth, we will consider the case when the initial distribution is lognormal, because this situation gives full meaning to the model in relation to the comments made above.

In this case (see Appendix B for a more detailed development) by maximizing the likelihood function it can be derived that the ML estimates of μ 1 and σ 2 (ln x i1μ 1 ) 2 , while the ML estimates for D = 1/c and A = e -k result of the solution of the system of equations (B.11) for the particular case of equally spaced data, that is t ijt i,j-1 = h, i = 1, . . . , d; j = 2, . . . , n i . This system has no explicit solution and must be dealt with by numerical methods. Once these estimates are obtained, the corresponding for b and σ 2 are b = b A, D and σ 2 = σ 2 A, D from (B.9). Finally, the ML estimate of any parametric function expressed in terms of b, c, k and σ 2 , as for example the mean function, is calculated by applying Zehna's theorem.

Numerical aspects

The system of equations (B.11) is quite complex, which makes its resolution difficult, especially when the sample is large. For this reason, it is necessary to make use of numerical procedures, most of which need an initial solution to be applied.

In order to obtain a good initial solution, we propose some alternatives based on the sample information provided by the observed trajectories of the process. In this line, we will distinguish several cases depending on the parameter b and will use some expressions derived from the curve.

Firstly, and taking into account that A = e -k and D = 1/c, from (4) we have

x 0 B ∞ = 1 - A t 0 D b , from where
D = A t 0 1 - x 0 B ∞ 1/b • (14)

A c c e p t e d m a n u s c r i p t

On the other hand, the inflection occurs at

t I = ln(bc) k = ln(D/b) ln A , so A = b 1 - x 0 B ∞ 1/b 1/(t 0 -t I ) . (15)
Now, we will distinguish the cases in which b is known and unknown. In each case we will consider two situations according to the inflection time being visualized or not.

• b known.

Firstly, let us suppose that t I > t 0 (i.e. b > e k t 0 /c). In such a case, and since 1 -

1 b b
is the quotient between the value of the curve at t I and the upper bound, we can approximate t I , for each trajectory, taking the first time instant at which the sample path exceeds

x i,∞ 1 -1 b b
, where x i,∞ is the upper bound for the i-th trajectory (i = 1, . . . , d), and, finally, consider the mean of these values. Obviously, this procedure can be applied if the limit value is known, at least in a approximate form. In this sense, usually the last value x i,n i of each trajectory is taken as x i,∞ .

Once t I is approached, the initial values for A and D are obtained from ( 14) and [START_REF] Rafail | A simple and precise method for fitting a Von Bertalannfy growth curve[END_REF], respectively. To this end, the value x 0 /B ∞ is approached by considering the mean of the values x i1 /x i,n i , i = 1, . . . , d. If b = 1 or ln(c)/k < t I ≤ t 0 (i.e., 1 < b ≤ e k t 0 /c), the inflection time cannot be guessed from the observation of the sample paths. Nevertheless, expression [START_REF] Lv | Stochastic Von Bertalanffy models, with applications to fish recruitment[END_REF] provides a relationship between A and D, or equivalently between k and c, concretely c = αe kt 0 where α = 1 -(x 0 /B ∞ ) 1/b . From this expression the curve remains

B(t) = B ∞ 1 -αe -k(t-t 0 ) b
and presents only an unknown parameter, k. We propose, for each sample path, to calculate k from a least square fit to the previous curve. For this, the value x 0 /B ∞ is approached in the same way described in the previous case, whereas B ∞ is approached by means of the values x i,n i , i = 1, . . . , d. Finally, the initial value for A is the exponential of the mean of these estimations.

• b unknown.

If for each sample path the inflection time can be visualized, or guessed, we propose to find an initial value for b as follows: firstly, we calculate an approximate value for the time at which the inflection occurs, for example by examining the sample paths, and taking the mean of the values

A c c e p t e d m a n u s c r i p t

x i,t I i /x i,∞ , namely x I , where t I i and x i,∞ are, respectively, the inflection time and the upper bound for the i-th sample path (approached by the last value of the trajectory). Once this value is obtained, the initial one for b is the solution of the equation

x I = 1 - 1 b b .
After calculating the initial value for b, the corresponding for A and D are obtained from ( 14) and ( 15) respectively, taking into account the same remarks made before.

If the inflection point cannot be visualized from the sample paths, that is ln(c)/k ≤ t I ≤ t 0 , we propose a slight modification of the previous procedure which is based on the increasing behavior of the von Bertalanffy curve. In fact, in such a case, and for b > 1, the following relation holds Once the value of b is obtained, we proceed as in the case when b is known.

1 - 1 b b = f (t I ) B ∞ ≤ x 0 B ∞ • Thus,

Simulation study

In this section we present several examples in order to validate the estimation procedure previously developed, together with the strategy showed for establishing the initial solution of the system of equations that must be solved. To this end, we have simulated sample paths from the von Bertalanffy diffusion process following the algorithm derived from the numerical solution of the stochastic differential equation associated with the process (see Kloeden et al. [START_REF] Kloeden | Numerical solution of SDE through computer experiments[END_REF]). In our case, the algorithm is A c c e p t e d m a n u s c r i p t

x n+1 = x n ⎡ ⎣ 1 + h bck e k tn -c - σ 2 2 + h 2 2 ⎛ ⎝ (bck 2 ) bc -e k tn (e k tn -c) 2 - σ 2 bck e k tn -c + σ 4 4 ⎞ ⎠ + σ 1 + hbck e k tn -c - hσ 2 2 1 + σ 2 Z 1n Z 1n + σ 3 6 Z 3 1n + σ 4 24 Z 4 1n ,
for n = 0, . . . , N -1, where N is the number of simulated values, h is the integration step, x n = X(t n ) and Z 1n is a normal variable of zero mean and variance h.

We have considered four cases in which d = 100 sample paths have been simulated. Each trajectory has been simulated with 301 data starting at t 0 = 0, taking h = 0.1 as the integration step and an initial lognormal distribution Λ(3, 0.2). Table 1 shows the values of the others parameters of the process, together with the theoretical inflection time instant. Note that in only in the first example the inflection time instant can be guessed. In order to make the subsequent inference we have considered, in each case, 31 data with t ijt i,j-1 = 1, i = 1, . . . , 100, j = 1, . . . , 31. The graphics in Figure 1 show the mean of the simulated sample paths for each case.

As concerns the estimation procedure, each of the examples has been treated considering the case b known and unknown.

• For example 1 (in which the inflection occurs inside the considered time interval), if b is known, the procedure aforementioned leads to 2.5 as approximate value or the inflection time (the true value is 2.35). We must note that the time instants considered are discrete; for this reason, when a time instant is found in which the value

x i,∞ 1 -1 b b
is exceeded, we have considered as the inflection time the mean between this time and the previous. Finally, from ( 14) and [START_REF] Rafail | A simple and precise method for fitting a Von Bertalannfy growth curve[END_REF] we have obtained the initial values for A and D and, therefore, those for c and k. In addition, the value x 0 /B ∞ has been approximated by 0.04042, whereas the true value is 0.04.

In the case b unknown, we have considered t I = 3 as the inflection time. From this value, and following the procedure described above, we have obtained b 0 = 1.22789, and then the others initial values as we have previously remarked. Table 2 summarizes the results obtained for this example. In the last case, the choice of t I is arbitrary and it is based on the visualization of the sample paths. For this reason we have considered other possible values for t I and compared the results with the previous one. Specifically, the values t I = 2, 4 have been also considered without changes in the estimated values (see table 3). • For example 2 we must remark that the inflection cannot be guessed (in fact, there is not inflection in the model since b = 1). In this example, the initial solution coincides in both cases, that is, when b is known and unknown because of the value obtained for x b is 0.40121 > e -1 . Therefore, when b is unknown, the procedure described leads to consider b = 1, from which follows that the initial solution for k (and then for A) and D is the same. Obviously, the final estimation of the parameters does not coincide in both cases since the likelihood equations differ. Table 4 contain the results for this example. exists but can not be seen. In this case we must remark that when b is unknown, its initial value is calculated by solving the equation

x b = 1 -1 b b since x b is now 0.20807 < e -1
. Thus, this equation has only one solution, b 1 = 1.59925, which is an upper bound for b. Later, we have taken a random value in the interval (1, b 1 ) as the initial value for b (Table 5 includes the results). Finally, the procedure for taking the initial value for b leads to consider a sensibility analysis of the results. To this end, we have considered other random values in the interval in order to validate the procedure followed (see table 6). However, no differences have been found, so the estimation of the parameters is robust to the choice of b. • The last example shows another situation in which the inflection exists but occurs before t 0 , so it can not be seen. The difference between this example and the previous is that, in this case, when b is considered unknown x b = 0.51208 > e -1 and then we take b 0 = 1 as the initial value for b. Please observe that, unlike in example 2, the initial values are different according to b being known or unknown. The results are included in Table 7.

A c c e p t e d m a n u s c r i p t

A c c e p t e d m a n u s c r i p t 6 Application to real data

The following application is based on the studies developed by J. Chong and M. Aguayo [START_REF] Chong | Edad y crecimiento del pez espada (Xiphias gladius Linnaeus, 1758) en el Pacífico suroriental (Diciembre 1994-Septiembre 1996)[END_REF] on some aspects related to the swordfish age and growth in the Southeastern Pacific. Specifically, we have considered data about the mean weight to age in both sexes, as shown in Figure 2. Table 8 shows the estimated values of the parameters of the model, together with the initial solution obtained by means of the strategy proposed before. To illustrate the predictive capability of the model, maximum likelihood estimates for the parameters using data for the age 1 to 16 were calculated, and we predicted the value for 17. The estimates of the parameters, and the initial solution, are summarized in Table 9, whereas Table 10 shows the observed values at age 17, point prediction given by the estimation of the mean and mode of X( 17)|X(16) = x 16 , as well as an interval prediction provided by the estimation of 0.025 and 0.975 quantiles of the variable above. One advantage of using growth stochastic models, like the one shown here, is the possibility of studying some questions related to their evolution, such

A c c e p t e d m a n u s c r i p t

as the achievement of a certain value or a change in a behavior pattern. For example, the time an animal takes to reach the minimum size at which it can be sold for consumption, or the time a population takes to reach a certain size, can be of great interest and can be formulated as a problem of first-passagetime through a constant boundary S, that is the time variable

T S,x 1 =Inf t≥t 1 {t : X(t) > S|X(t 1 ) = x 1 }
for which its probability density function must be calculated. In order to illustrate this question, we have considered two boundaries, S = 300 and S = 350. For the first case the value S is achieved inside the time interval considered, whereas the second will be attained later (the estimated upper bound for the model is 405.755).

As regards the obtaining of the probability density function, it is the solution of a Volterra integral equation of the second kind that must be solved by numerical procedures since a closed-form solution is not available in this case (see [START_REF] Gutiérrez | First-passagetime densities for time-non-homogeneous diffusion processes[END_REF]). To this end we have employed the methodology developed in [START_REF] Román | First-passage-time location function: Application to determine first-passage-time densities in diffusion processes[END_REF].

Conclusions

A new stochastic diffusion model has been presented for the modeling of animal growth with a behavior that fits a von Bertalanffy growth pattern. The main innovation from other existing models (see [START_REF] Gudmundsson | Stochastic growth[END_REF], [START_REF] Lv | Stochastic Von Bertalanffy models, with applications to fish recruitment[END_REF], [START_REF] Russo | Lévy processes and stochastic von Bertalanffy models of growth, with appliccation to fish population analysis[END_REF]) is that it allows for modeling both length and weight for some animal species in isometric and The new model is a diffusion process determined by the solution of a differential stochastic equation, which can be obtained by including a noise term in the ordinary differential equation associated with the respective deterministic model. Thus, and in the absence of noise, the growth variable increases according to [START_REF] Capocelli | Growth with regulation in random environment[END_REF] and the sample paths of the process show such behavior patterns (according to a generalized von Bertalanffy curve), save for the presence of random disturbances (whose magnitude depends on that of the noise). In addition, the fact that the mean function of the process is also a type (2) curve, makes the model particularly interesting for real data applications with fitting and forecasting purposes.

In addition to the mean function, other parametric functions prove useful for fitting and forecasting aims. In the case of point fitting and prediction, the mode function provides, for each time, the most probable value for the growth variable. For interval predictions, through consideration of the percentile function, intervals containing the growth variable for each time (with a specific probability), can be obtained.

Once the diffusion process has been defined, the maximum likelihood estimation of its parameters is done (and, consequently, that of the parametrical functions of interest). However, the likelihood equations obtained do not have an explicit solution and classical numerical procedures are required. The convergence of such procedures may depend on the choice of good initial solutions.

For this reason, one of the contributions of the present paper is a proposal for a strategy for the search of such initial solutions in several potential situations. The cases of b known and unknown have been considered and also, in order to use the information on the characteristics of the model provided by the sample data, we have distinguished two possibilities regarding the visualization of the inflection point in the sample paths. Such strategy is validated by means of simulated examples that consider all the situations described.

Finally, one of the main advantages of the use of stochastic models is, in addition to its fitting and forecasting capabilities, the chance to study time variables that affect the growth process. In an application to real data of mean weight of swordfish, the capability of the model for fitting and predicting is shown, as well as for studying time variables of interest like the first-passagetime through constant boundaries.

A c c e p t e d m a n u s c r i p t

1 τ E X (n+1)τ -X nτ |X nτ = x = bck e knτ -c x, 1 τ E X (n+1)τ -X nτ 2 |X nτ = x = σ 2 x 2 , 1 τ E X (n+1)τ -X nτ 2+p |X nτ = x = x 2 τ o(τ ), ∀p ∈ N, (A.2)
in such a way that when τ approaches zero and n → ∞, when the condition nτ = t, then X nτ converges to the diffusion process with infinitesimal moments provided by [START_REF] Gudmundsson | Stochastic growth[END_REF].

B Appendix B

The transition p.d.f. ( 12) can be rewritten as

f (x ij , t ij | x i,j-1 , t i,j-1 ) = 1 x ij 2πσ 2 (t ij -t i,j-1 ) × exp ⎛ ⎜ ⎝- 1 2 ln x ij x i,j-1 -b ln D-A t ij D-A t i,j-1 + σ 2 2 (t ij -t i,j-1 ) 2 (t ij -t i,j -
L x ij ( μ 1 , σ 2 1 , A, D, b, σ 2 ) = - n 2 ln(2π) - d 2 ln σ 2 1 - n -d 2 ln σ 2 - d i=1 ln x i1 - 1 2σ 2 1 d i=1 [ln x i1 -μ 1 ] 2 - d i=1 n i j=2 ln x ij - 1 2 d i=1 n i j=2 ln(t ij -t i,j-1 ) - 1 2 d i=1 n i j=2 [ln x ij x i,j-1 -b ln D-A t ij D-A t i,j-1 + σ 2 2 (t ij -t i,j-1 )] 2 (t ij -t i,j-1 )σ 2 ,
from which the ML estimates of μ 1 and σ x i,j-1b ln D-A t ij D-A t i,j-1 + σ 2 2 (t ijt i,j-1 ) t ijt i,j-1 × t i,j-1 A t i,j-1 -1 (D -A t ij )t i,j A t i,j -1 (D -A t i,j-1 ) (D -A t ij ) (D -A t i,j-1 ) = 0, (B.1)

d i=1 n i j=2 ln x ij
x i,j-1b ln D-A t ij D-A t i,j-1 + σ 2 2 (t ijt i,j-1 ) t ijt i,j-1 × ln D -A t ij D -A t i,j-1 = 0, (B.2) x i,j-1 ln D-A t ij D-A t i,j-1 t ijt i,j-1 = 0.(B.4)

d i=1 n i j=2 ln x ij x i,j-1 -b ln D-A t ij D-A t i,j-1 + σ 2 2 (t ij -t i,j-1 ) t ij -t i,j-1 × A t ij -A t i,j-1 (D -A t ij ) (D -A t i,j-1 ) = 0, (B.
This system of equations does not have an explicit solution, so numerical procedures are required to find it. In the particular case t ijt i,j-1 = h, i = 1, . . . , d; j = 2, . . . , n i , equations (B.1) to (B.4) remain

D(1 -A h ) 2X A,D 2

A c c e p t e d m a n u s c r i p t

  at successive time instants can be time-consuming and expensive. For this reason, it is common to consider a single sample path composed of the mean of the recapture data in each time instant.

- 1 ,

 1 and in an approximate form, we can consider that 1 -1 b b < x b , where x b = Mean i=1,...,d x i,1 /x i,n i . Furthermore, 1 -∀b ≥ 1, from which we can establish the following strategy in order to take an initial value for b in this situation: If x b < e -1 , then the solution of the equation x b = 1 -1 b b provides an upper bound for b, namely b 1 . In such a case, we propose to consider an initial value for b randomly chosen from an uniform distribution in the interval (1, b 1 ). If x b ≥ e -1 , then the previous equation has no solution, so we propose to consider b = 1 as initial value.

Figure 1 .

 1 Figure 1. Mean of the simulated sample paths

•

  As example 3 concerns, we have another situation in which the inflection

Figure 2 .

 2 Figure 2. Mean weight to age for swordfish in both sexes.

A c c e p t e d m a n u s c r i p tFigure 3

 3 Figure 3 shows the observed values, the estimation of the E[X(t)|X(t 1 ) = x 1 ] function and the estimation of the E[X(t i )|X(t i-1 ) = x i-1 ] values (joined by a solid line for better visualization), i = 1, . . . , 17.

Figure 3 .

 3 Figure 3. Observed values and estimations of the trend function and the E[X(t i )|X(t i-1 ) = x i-1 ] values, i = 1, . . . , 17.

Figure ( 4 )

 4 shows the probability density functions for the two cases considered for which, moreover, we can point out that the means of the variables firstpassage-time through the constant S are 14.138 and 18.903, respectively.

Figure 4 .

 4 Figure 4. First-passage probability densities for the weight of the swordfish through a constant boundary: (a) S = 300; (b) S = 350.

A c c e p t e d m a n u s c r i p t

  allometric situations. This is done by considering a generalized von Bertalanffy growth curve including a parameter b ≥ 1, and whose particular case with b = 1 is the one considered in previous stochastic models.

1 )σ 2 ⎞

 12 D = 1/c and A = e -k . Denoting n = d i=1 n i , the log-likelihood function of the sample is ln

(ln x i1 -μ 1 ) 2 ,

 2 whereas the corresponding to A, b, D and σ 2 follow from the solution of the system of equations

  t i,j-1 ) + 4σ 2 (nd) -

2 -2bW A,D 3 + σ 2 hσ 4 3 + 8bY A,D 2 - 1 . 2 Y A,D 1 -X A,D 1 Y A,D 2 X A,D 3 Y A,D 1 -X A,D 1 Y 2 Y A,D 3 -X A,D 3 Y A,D 2 X A,D 3 Y A,D 1 -X A,D 1 Y A,D 3 . (B. 10 )D( 1 - 2 - 3 + 2 Y A,D 3 + 2 -A c c e p t e d m a n u s c r i p t

 232432121123112332311310123232 , * -2bX A,D 3, * + σ 2 h X A,D 1, * +hA h-1 2W A,D h 2 (nd) + 4σ 2 h(nd) -4b 2 Y A,D 4Z = 0, (B.8)where, and denotingS A,D ij = (D -A t ij ) (D -A t i,j-1 ) and T A,D ij = ln D-A t ij D-A t i,j-1 ,After some algebra, from equations (B.6) and (B.7), we obtainb A,D = X A,DReplacing these expressions in (B.5) and (B.8) the following system of equations appearsA h ) X A,D 2, *b A,D X A,D 3, * + C A,D X A,D 1, * +hA h-1 W A,D b A,D W A,D d) C A,D 2 + 2(nd)C A,Db A,D 2b A,D Y A,D Z = 0.As we pointed out in the introduction, in some cases the parameter b is known. For example, in the study of the length of fish b is equal to one, whereas in the case of weight in isometric growth b = 3. In such a case, that is, when b

Table 1

 1 Parameters chosen for the simulation of the sample paths

	Example b	c	k	σ	t I
	1	2 0.8 0.2 0.01 2.35
	2	1 0.6 0.5 0.01 -1.02
	3	1.3 0.7 0.4 0.01 -0.23
	4	3 0.2 0.3 0.01 -1.70

Table 2

 2 Initial and estimated values for example 1

				b known		b unknown
	Parameter True value Initial solution Estimated value Initial solution Estimated value
	b	2					1 .22789	1.97997
	c	0.8		0 .79894	0.79994		0.92668	0.80339
	k	0.2		0 .18747	0.20024		0.04305	0.19847
	σ	0.01			0.00991			0.00991
	Table 3
	Sensibility analysis for example 1 (b unknown).
	t I	b 0	c 0	k 0	b	c	k	σ
	2 1.09777 0.94621 0.01899 1.97997 0.80339 0.19847 0.00991
	4 1.49086 0.88375 0.06894 1.97997 0.80339 0.19847 0.00991
	Table 4
	Initial and estimated values for example 2
				b known		b unknown
	Parameter True value Initial solution Estimated value Initial solution Estimated value
	b	1					1	1 .00201
	c	0.6		0 .59878	0.60125		0.59878	0.60428
	k	0.5		0 .50850	0.49515		0.50850	0.49354
	σ	0.01			0.00996			0.01011

Table 5

 5 Initial and estimated values for example 3

				b known		b unknown
	Parameter True value Initial solution Estimated value Initial solution Estimated value
	b	1.3				1 .37891	1.30028
	c	0.7	0 .70107	0.70075		0.67968	0.70067
	k	0.4	0 .39513	0.39842		0.40180	0.39847
	σ	0.01		0.01017			0.01017
	Table 6
	Sensibility analysis for example 3 (b unknown).
	b 0	c 0	k 0	b	c	k	σ
	1.11985 0.75385 0.85843 1.30028 0.70067 0.39847 0.01017
	1.17977 0.73568 0.73702 1.30028 0.70067 0.39847 0.01017
	1.2397 0.71813 0.62710 1.30028 0.70067 0.39847 0.01017
	1.29962 0.70118 0.52597 1.30028 0.70067 0.39847 0.01017
	1.35955 0.68484 0.43135 1.30028 0.70067 0.39847 0.01017
	1.41947 0.66909 0.34097 1.30028 0.70067 0.39847 0.01017
	1.4794 0.65393 0.25291 1.30028 0.70067 0.39847 0.01017
	1.53932 0.63934 0.17442 1.30028 0.70067 0.39847 0.01017

Table 7

 7 Initial and estimated values for example 4

				b known	b unknown
	Parameter True value Initial solution Estimated value Initial solution Estimated value
	b	3			1	3 .10916
	c	0.2	0 .19995	0.19998	0.40073	0.19361
	k	0.3	0 .29958	0.30198	0.26639	0.30305
	σ	0.01		0.00980		0.00980

Table 8

 8 Estimation of the parameters of the model using all swordfish data. ) = x 1 ], where t 1 is the first observation time instant and x 1 is the initial value of the sample path. Nevertheless, and with the aim of obtaining a best fit to the observed data, in practice it is usual to consider the values E[X(t i )|X(t i-1 ) = x i-1 ], i = 1, . . . , 17, where x i-1 is the observed value at t i-1 .

	Parameter Initial solution Estimated value
	b	4.2379	3.0631
	c	0.5901	0.7305
	k	0.1833	0.1451
	σ		0.0003
	Since the model considered verifies that E[X(t)|X(t 0 ) = x 0 ] is a generalized
	von Bertalanffy curve, it is obvious that a good fit to the data will be provided
	by the function E[X(t)|X(t 1

Table 9

 9 Estimation of the parameters of the model without the last data.

		Parameter Initial solution Estimated value
		b	5.2573	3.0589
		c	0.5067	0.7310
		k	0.1959	0.1450
		σ		0.0003
	Table 10
	Observed and predicted values for the mean weight of swordfish at 17-th year.
	Observed	Estimated	Estimated	Estimated
	value	conditional mean conditional mode conditional quantiles
	333.2	333.331	333.445	(333.092,333.571)

A c c e p t e d m a n u s c r i p t 8 Acknowledgements

The authors are very grateful to the anonymous referees whose comments and suggestions have contributed to improve this paper. This work was supported in part by the Ministerio de Educación y Ciencia, Spain, under Grants MTM2008-05785/MTM and HI2007-0034, and by the Consejería de Innovación, Ciencia y Empresa (Junta de Andalucía), Spain, under Grant P06-FQM-02271.

A Appendix A

Another way to develop the diffusion process proposed in this paper is in the line shown by Albano and Giorno [START_REF] Albano | A stochastic model in tumor growth[END_REF]. The procedure is based on the discrete version of (8),

that approaches [START_REF] Gutiérrez | Approximate and generalized confidence bands for some parametric functions of the lognormal diffusion process with exogenous factors[END_REF] when n → ∞ and τ → 0, with nτ = t.

Now we introduce random environment as follows: the relative intrinsic change bckτ / e knτc in the interval [nτ, (n + 1)τ ), n = 0, 1, . . . can bee seen as the mean value of a sequence of independent (not identically distributed) Bernoulli variables, Z nτ , n = 0, 1, . . ., verifying

where σ > 0 is a constant that evaluates the width of the environmental fluctuations. Furthermore, the moments of Z nτ are

In this way, we randomize the model by replacing in the term bckτ / e knτc by Z nτ (A.1), thus obtaining