## Abstract

By using functional integral methods, we determine a computable evolution equation for the joint response-excitation probability density function of a stochastic dynamical system driven by coloured noise. This equation can be represented in terms of a superimposition of differential constraints, i.e. partial differential equations involving unusual limit partial derivatives, the first one of which was originally proposed by Sapsis & Athanassoulis. A connection with the classical response approach is established in the general case of random noise with arbitrary correlation time, yielding a fully consistent new theory for non-Markovian systems. We also address the question of computability of the joint response-excitation probability density function as a solution to a boundary value problem involving only one differential constraint. By means of a simple analytical example, it is shown that, in general, such a problem is undetermined, in the sense that it admits an infinite number of solutions. This issue can be overcome by completing the system with additional relations yielding a closure problem, which is similar to the one arising in the standard response theory. Numerical verification of the equations for the joint response-excitation density is obtained for a tumour cell growth model under immune response.

## 1. Introduction

The relevance of coloured noise in physical systems has been recognized in many different applications, including stochastic resonance in sensory neurons (Nozaki *et al.* 1999), chemical excitable systems (Zhong & Xin 2001; Beato *et al.* 2008), structural dynamics (Bolotin 1969; Shlesinger & Swean 1998; Li & Chen 2009), tumoural cell growth (Fiasconaro *et al.* 2006; Wang 2009; Zeng & Wang 2010), spin relaxation in magnetic phenomena (Atxitia *et al.* 2009), statistical properties of dye lasers (Luo *et al.* 2001; Jin *et al.* 2005) and optical instabilities (Jung & Hänggi 1988; Zhou 2009). Owing to the inherent complexity of these systems, a comprehensive mathematical description is not viable in practice, and therefore their study often relies on reduced-order models consisting of nonlinear ordinary differential equations subject to uncertainty in physical parameters, initial conditions or including random forcing terms with non-trivial correlation structure. The statistical properties associated with the solution to these initial value problems can be computed according to different stochastic methods. Well-known approaches are generalized polynomial chaos (Ghanem & Spanos 1998; Xiu & Karniadakis 2002, 2003), multi-element generalized polynomial chaos (Wan & Karniadakis 2006; Venturi *et al.* 2010), multi-element and sparse grid adaptive probabilistic collocation (Foo & Karniadakis 2008, 2010), high-dimensional model representations (Rabitz *et al.* 1999; Li *et al.* 2002; Yang *et al.* submitted), stochastic biorthogonal expansions (Venturi 2006, 2011; Venturi *et al.* 2008) and proper generalized decompositions (Chinesta *et al.* 2010; Nouy 2010*b*).

An interesting new approach for the complete probabilistic description of a stochastic system driven by coloured noise was recently proposed by Sapsis & Athanassoulis (2008). The key idea—inspired by the work of Beran (1968)—was to perform an analysis on the extended probability space consisting of the joint *response-excitation* statistics. In particular, the Hopf equation (Lewis & Kraichnan 1962; Rosen 1971; Klyatskin 1974) governing the dynamics of the joint response-excitation characteristic functional of the system was reduced to a *differential constraint* involving the one-point response-excitation probability density function. This differential constraint was then supplemented with additional marginal compatibility conditions and other constitutive relations in order to obtain a closed system of equations. This led to the formulation of a new theory, which was shown to be consistent with standard stochastic approaches such as moment equations or Fokker–Planck equations. In this way, the closure problem arising in the standard response approach for non-Markovian systems (Stratonovich 1967; Moss & McClintock 1995) was apparently overcome. However, the transport equation obtained for the joint response-excitation density (Sapsis & Athanassoulis 2008, eq. (5.7a)) involves an unusual partial differential operator, i.e. a limit partial derivative, which makes the question of well-posedness and computability not clear. The purpose of this paper is to study in more detail this fundamental question and therefore assess the effectiveness of the response-excitation theory for the simulation of stochastic systems driven by coloured noise. To this end, we will employ a functional integral method that simplifies considerably the derivation of the differential constraints given by Sapsis & Athanassoulis (2008), and it also allows to generalize them for systems having smooth nonlinearities of non-polynomial type.

This paper is organized as follows. In §2, we determine the set of equations arising from the response-excitation theory for a simple first-order nonlinear ordinary differential equation driven by purely additive random noise. We also discuss the question of computability and well-posedness of the boundary value problem for the joint response-excitation density of the system. The connection between the response-excitation approach and the classical response approach is established in §3 for random noise with arbitrary correlation time. In §4, we review useful formulae to separate the correlation structure between two functionals of the random forcing process. These formulae will allow us to deal effectively with weakly coloured noise approximations of stochastic systems driven by Gaussian forces. Ordinary differential equations involving random parameters are treated in §5. In §6, we present a numerical application to the tumour cell growth model recently proposed by Zeng & Wang (2010). Finally, the main findings and their implications are summarized in §7. We also include two brief appendices dealing, respectively, with the formulation of differential constraints for the joint response-excitation probability density of a nonlinear pendulum driven by a random torque (appendix A) and with the theory of small correlation time approximations (appendix B).

## 2. An evolution equation for the joint response-excitation probability density

In this section, we develop a systematic methodology to determine an equation for the joint response-excitation probability density function of a stochastic dynamical system driven by coloured noise. To this end, let us first consider the following simple prototype problem:
2.1
where *f*(*t*;*ω*) is a smooth coloured random noise while *g*(*x*,*t*) is a nonlinear function, which is assumed to be Lipschitz continuous in *x* and continuous in *t*.^{1} The solution to (2.1) (when it exists) is a *non-local* functional of the forcing process *f*, which will be denoted as *x*_{t}[*f*]. Similarly, we use the shorthand notation *f*_{s} to identify the random variable *f*(*s*;*ω*), i.e. the stochastic process *f* at time *s*. We also assume that the problem (2.1) admits the existence of the joint response-excitation probability density function, i.e. the joint probability density function of the *response process* *x*_{t}[*f*] at time *t* and the *excitation process* *f*_{s} at time *s*. Such a probability density has the following functional integral representation:
2.2
where *Q*[*f*] denotes the probability density functional of the random forcing (Fox 1986, p. 467) while is the functional integral measure (Phythian 1977; Jouvet & Phythian 1979; Jensen 1981). The representation (2.2) is based on the following finite dimensional result (Khuri 2004, eq. (15))
2.3
where *Q*(*f*(*t*_{1}),…,*f*(*t*_{N})) denotes the joint probability density of the forcing process at times *t*_{1},…,*t*_{N} (i.e. a random vector), while *x*_{ti}(*f*(*t*_{1}),…,*f*(*t*_{N})) is a nonlinear mapping from into representing the response process at time *t*_{i} as a function of the forcing process at times *t*_{1},…,*t*_{N}. The functional integral (2.2) is defined as the limit of equation (2.3) as *N* goes to infinity. In this sense, in equation (2.2) we recognize the standard definitions of *Q*[*f*] and as given, e.g. by Phythian (1977) and Langouche *et al.* (1979)
2.4
2.5
2.6
for an arbitrary discretization of the integration period into *N* points.

An evolution equation for the joint response-excitation probability density of the system (2.1) can be determined by differentiating equation (2.2) with respect to *t*. This yields
2.7
Taking the limit for gives the following result, first obtained by Sapsis & Athanassoulis (2008) using a Hopf characteristic functional approach
2.8
This equation looks like a closed evolution equation for the joint response-excitation probability density function of the stochastic dynamical system along the direction *s*=*t*. In the work of Sapsis & Athanassoulis (2008), equation (2.8) was also accompanied with an initial condition and with the *marginal compatibility condition*:^{2}
2.9
expressing the fact that the evolution of the joint response-excitation density has to be consistent with the evolution of the excitation density . Note that since the stochastic process *f*(*t*;*ω*) is given, is a known function. In addition, the joint density has to satisfy the following two obvious, yet essential, *constitutive conditions*^{2}
2.10
The system (2.8)–(2.10) was proposed as a computable set of equations describing the evolution of the joint response-excitation probability density function of dynamical systems driven by coloured noise. Note that in the extended probability space consisting of the joint response-excitation statistics, the standard closure problem arising, e.g. in the classical coloured noise master equation (Hänggi 1985) seems to be overcome. However, the presence of the limit partial derivative at the left-hand side of equation (2.8) should warn us on the fact that we are not dealing with a standard partial differential equation. Indeed equation (2.8) is rather a *differential constraint* (Venturi & Karniadakis 2011) to be satisfied by the joint response-excitation probability density function of any solution to equation (2.1) along the line *s*=*t*. Preliminary insight into the question of well-posedness of the boundary value problem (2.8)–(2.10) can be gained by expanding equation (2.7) for *s* in the neighbourhood of *t*. This gives us a first-order correction to the differential constraint (2.8), which ultimately leads to a standard partial differential equation for the joint density on the infinitesimal strip
2.11
To this end, let us assume that the random forcing process is differentiable in time and expand it into a Taylor series around *s* as
2.12
A substitution of this expansion back into equation (2.7) yields, after simple mathematical manipulations,
2.13
This equation holds for all *t*≥*t*_{0} and for , i.e. *s* in the neighbourhood of *t* (see (2.11)). Analysis of equation (2.13) clearly shows that the derivatives and are *coupled*. In other words, in order to compute the joint response-excitation probability density function in the neighbourhood of *s*=*t*, we need an additional expression for . This is also the case when we take the limit . In fact, by using the obvious identity
2.14
we can see that the dynamics of the joint response-excitation probability density function in the direction *s*=*t* (i.e. a directional derivative) can be represented as a superimposition of two differential constraints: the first one is equation (2.8), while the second one is
2.15
This yields the partial differential equation
2.16
which is the complete evolution equation governing the dynamics of the joint response-excitation probability density function of the system. Equation (2.16) can be evaluated further if one has available an expression for the average or, equivalently, 〈*δ*(*a*−*x*_{t}[*f*])*δ*(*b*−*f*_{s})〉. Such expression involves non-local functionals of the random forcing process *f* and it will be provided in §4. As we shall see in the next subsection, if we do not include the limit derivative (2.15)—i.e. the differential constraint complementary to equation (2.8)—within the set of equations, then the system (2.8)–(2.10) turns out to be *undetermined*, in the sense that it possibly admits an infinite number of solutions.^{3}

### (a) Example: ill-posedness of the boundary value problem (2.8)–(2.10)

Let us consider the following trivial first-order dynamical system driven by a smooth random force *f*(*t*;*ω*):
2.17
For the purposes of the present example, it is enough to set , where *ξ*(*ω*) is a Gaussian random variable. Let us also assume that the initial state of the system *x*_{0}(*ω*) is Gaussian distributed with zero-mean, and that *x*_{0}(*ω*) is independent of *ξ*(*ω*). The analytical solution to (2.17) is obviously
2.18
This allows us to obtain the joint probability of *x*(*t*;*ω*) and *f*(*s*;*ω*) by using the classical mapping approach (Papoulis 1991, p. 142). Specifically, we consider the following mapping between the random variables (*ξ*(*ω*),*x*_{0}(*ω*)) and (*x*(*t*;*ω*),*f*(*s*;*ω*)):
2.19
where
2.20
This yields the joint response-excitation probability density function
2.21
It is straightforward to verify that (2.21) satisfies
2.22
which is the equation arising from equation (2.8) by setting *g*(*a*,*t*)=−*a*. Also, the marginal compatibility condition (2.9) as well as the constitutive relations (2.10) are obviously verified. However, if we set *C*(*t*)=*α*(*t*) in (2.21), where *α*(*t*) is an *arbitrary function* such that *α*(*t*_{0})=0, then we easily see that we still have a probability density function satisfying equations (2.22), (2.9) and (2.10) and the initial condition. This suggest that this boundary value problem is not well-posed, in the sense that it admits an infinite number of solutions. In addition, we remark that setting *C*(*t*)=*α*(*t*) in (2.21) is not the only degree of freedom we have, as other solutions can be constructed. These observations provide a definitive answer to the question raised by Sapsis & Athanassoulis (2008, p. 295), regarding the solvability theory of the system (2.8)–(2.10). The multiplicity of solutions admitted by such system arises because the correlation structure between the response process *x*_{t} and the excitation process *f*_{s} was not properly taken into account in the formulation of the theory.

## 3. An evolution equation for the response probability density

In the past few decades, many researchers attempted to determine a closed equation describing the evolution of the response probability density of a stochastic system driven by coloured random noise (Stratonovich 1967; Fox 1986; Faetti & Grigolini 1987; Hänggi & Jung 1995; Moss & McClintock 1995). Perhaps, the first effective approach was developed by the school surrounding Stratonovich (1967). The starting point is the functional representation of the response probability density. For the system (2.1), we have
3.1
where the average is with respect to the joint probability functional of the excitation process and the initial state. Differentiation of (3.1) with respect to time yields
3.2
This equation can be evaluated further if one has available an expression for the average appearing in the last term at the right-hand side. Among different methods devised to represent such quantity, we recall small correlation time expansions (Stratonovich 1967; Dekker 1982; Fox 1986; Lindenberg & West 1983) (see also §4 and appendix B), cumulant resummation methods (Lindenberg & West 1983; Lindenberg *et al.* 1989), functional derivative techniques (Hänggi 1978*a*,*b*, 1985), path integral methods (Pesquera *et al.* 1983; Wio *et al.* 1989; McCane *et al.* 1990; Venkatesh & Patnaik 1993), decoupling approximations (Hänggi & Jung 1995) and operator projection methods (Grigolini 1981; Faetti & Grigolini 1987).

### (a) Consistency of the response-excitation theory with the classical response theory

It is important at this point to prove that the response-excitation theory is consistent with classical approaches for the response probability density function. This establishes a full correspondence between the two methods. To this end, let us first consider the differential constraint (2.8) and integrate it with respect to the variable *b* from to . This yields
3.3
Now, let us take a closer look at the last term at the right-hand side of equation (3.3). By definition (2.2), we have
3.4
If we substitute this result into equation (3.3), we obtain exactly equation (3.2). Therefore, we have shown that the differential constraint (2.8) is consistent with the classical response approach for random noise with arbitrary correlation time. This result extends the one obtained by Sapsis & Athanassoulis (2008) for systems driven by white noise.

Next, we consider the evolution equation (2.16). By following the same steps that led us to the consistency result just discussed, we can show that the classical response approach is also included in equation (2.16). In fact, if we perform an integration with respect to *b* from to , we see that the last term in equation (2.16) vanishes. This is due to the properties of the probability density function at . Therefore, both the differential constraint (2.8) and the full evolution equation (2.16) are consistent with the classical response theory.

## 4. Correlation between two functionals of the random forcing process

By using functional analysis, it is possible to obtain a representation of the correlation structure between two arbitrary functionals *F*_{t}[*f*] and *G*_{s}[*f*] of the random forcing process (where the subscripts in *F* and *G* denote the time at which those functionals are evaluated) in terms of a functional power series involving the cumulants of the forcing process itself. The general result was first determined by Klyatskin (1974) (see also Klyatskin 2005, p. 70) in an operator form as
4.1
where the symbol *δ*/*δ*(⋅) denotes a Volterra functional derivative (Volterra 1959; Beran 1968) and *Z*_{f} is the Hopf characteristic functional of the excitation process, i.e.
4.2
Later on, Bochkov & Dubkov (1974) obtained similar results, first for the correlation of two regular functionals of a Gaussian random process and subsequently for two regular functionals of two arbitrary processes (Bochkov *et al.* 1977). The general formula is
4.3
where denote the cumulants of *f*. We recall that equation (4.3) is valid for all *s*<*t*. Some care is needed if we consider functionals up to final time *s*=*t* (Hänggi 1978*a*; Hänggi & Jung 1995). This fact has been also pointed out by Klyatskin (2005, p. 71), who noticed that the expansion (4.3) does not always give the correct result in the limit . This means that the limiting procedure and the procedure of expansion in the functional power series sometimes do not commute.

In the particular case where *G*_{s}[*f*]=*f*_{s}, the cumbersome relation (4.3) simplifies to
4.4
This generalizes the well-known Furutsu–Novikov–Donsker formula
4.5
which is valid when *f* is a Gaussian process.^{4} A substitution of equation (4.4) with *F*_{t}[*f*]=*δ*(*a*−*x*_{t}[*f*]) into equation (3.2) yields the so-called *coloured noise master equation* pioneered by Hänggi (1978*a*,*b*). Such equation describes the evolution of the response probability of a first-order dynamical system subjected to additive zero-mean coloured random noise
4.6
The summation at the right-hand side of this equation involves the response process itself, i.e. *x*_{t}[*f*], through the response function of the system. Therefore, equation (4.6) is not in a closed form, meaning that it requires further manipulations in order to be computable. In the simpler case of zero-mean Gaussian forcing processes, equation (4.6) reduces to
4.7
where
4.8
Therefore, the problem of finding the evolution of the response probability in the presence of *Gaussian* forcing has been reduced to the evaluation of the *response function* *δx*_{t}[*f*]/*δf*_{s}, which is given as a functional Volterra derivative of the process *x*_{t}[*f*] with respect to *f*_{s}. A general framework for calculating such a response function is the operator approach of Martin *et al.* (1973) (see also Phythian 1977; Jouvet & Phythian 1979; Jensen 1981). However, for simple dynamical systems described by equation (2.1), the response function can be calculated directly from the equation of motion. In fact, if we functionally differentiate (2.1) with respect to *f*_{s}, we obtain
4.9
where we have denoted by . This is a first-order linear differential equation in *δx*_{t}[*f*]/*δf*_{s} whose general solution, corresponding to the initial condition *δx*_{t0}[*f*]/*δf*_{s}=0,^{5} is given by
4.10
where *Θ*(*t*−*s*) is a step function expressing causality (i.e. the force *f*_{s} acts on the path *x*_{t}[*f*] only when *s*≤*t*). At this point, we can write equation (4.7) as
4.11
This shows that the evolution of the response probability density theoretically involves the entire history of the response process *x*_{t}[*f*] from *t*_{0} to *t* in a functional integral form. Many approximate analytical methods have been developed in the past in order to obtain explicit closures of equation (4.6) or equation (4.11). In appendix B, we recall the small correlation time approximation method.

### (a) Separation of the joint response-excitation probability density function

The joint response-excitation probability density function (2.2) is defined as the average of the product of two functionals of the excitation process, namely *F*_{t}[*f*]=*δ*(*a*−*x*_{t}[*f*]) and *G*_{s}[*f*]=*δ*(*b*−*f*_{s}). We have already seen how to disentangle the correlation structure of two functionals of the random forcing process in terms of the statistical properties of the process itself (see equation (4.3)). In the particular case of Gaussian forcing processes, we have that *C*_{i}=0 for *i*>2, and therefore the general expression (4.3) simplifies to
4.12
This formula was first obtained by Bochkov & Dubkov (1974). By using equation (4.12), we can separate the joint response-excitation density (2.2). To this end, we simply substitute *F*_{t}[*f*]=*δ*(*a*−*x*_{t}[*f*]) and *G*_{s}[*f*]=*δ*(*b*−*f*_{s}) in (4.12) to obtain
4.13
where we have used the identities
4.14
Clearly, if *f* is white noise and *s*<*t*, then equation (4.13) implies the well-known result
4.15
This suggests that when we are dealing with a weak correlation structure in the noise, only the linear term in *C*_{2} should be retained in the expansion (4.13). In this case, the joint probability density function (2.2) can be approximated as
4.16
The first term in (4.16) represents the white noise approximation of when *s*<*t*. The second term can be evaluated explicitly by using either equation (B5) or equation (B7) obtained in appendix B. This yields a representation of the joint response-excitation probability density function, which is based on a time asymptotic closure of the correlation structure between the response and the excitation processes. For instance, in the case of weakly coloured exponential correlation function, we obtain
4.17
where ℓ and *D* are two positive parameters that control, respectively, the correlation time and the correlation amplitude of the noise (see appendix B). It can be shown that formula (4.17) satisfies the marginal compatibility condition (2.9).

### (b) Computable evolution equations for the response and the response-excitation probability density functions in a weakly coloured approximation

By using the small correlation time expansions obtained in appendix B, it is possible to obtain a closure approximation of equation (4.11) exploiting the fact that the response process, in this case, is nearly Markovian. To this end, we simply substitute, e.g. equation (B5) back into equation (4.11) to obtain
4.18
where ℓ and *D* are the correlation time and the correlation amplitude of *f* (see equation (B1)), respectively. Similarly, a substitution of equation (4.17) into the last term at the right-hand side of equation (2.16) results in the following integro-differential equation for the joint response-excitation probability density function of the system:
4.19
Let us recall that evolution equations based on small correlation time expansion are known to have some technical difficulties owing to possible negative diffusivity and absence of uniform convergence (Fox 1984; Fox 1986, p. 471; Hänggi 1989, p. 319).

## 5. Ordinary differential equations involving random parameters

The response-excitation theory developed in the previous sections can be applied also to equations involving random parameters. This class of problems obviously includes equations driven by random processes represented in terms of a Karhunen–Loève series. To fix ideas, let us consider the first-order system
5.1
where (*ξ*_{1}(*ω*),…,*ξ*_{M}(*ω*)) is a vector of random variables with known joint probability density function. For instance, the multivariable function *h* could be in the form
5.2
the second term being the finite-dimensional Karhunen–Loève representation of the random forcing appearing in equation (2.1). In order to determine an evolution of the joint response-excitation density associated with the system (5.1), it is convenient to identify as ‘excitation’ the random vector (*ξ*_{1}(*ω*),…,*ξ*_{M}(*ω*)). Thus, we look for an equation satisfied by the joint density
5.3
A time differentiation of (5.3) yields
5.4
5.5
i.e.
5.6
This is a possibly high-dimensional but closed and exact linear partial differential equation in two variables (*t* and *a*) and *M* parameters (*b*_{1},…,*b*_{M}) which can be solved numerically by exploiting recent advances in numerical methods for high-dimensional systems such as proper generalized decomposition (Chinesta *et al.* 2010; Nouy 2010*a*), sparse grid collocation (Foo & Karniadakis 2008, 2010) or functional ANOVA techniques (Rabitz *et al.* 1999; Cao *et al.* 2009; Yang *et al.* submitted). We remark that equations of the type (5.6) were first determined a long time ago by Dostupov & Pugachev (1957). More recently, Li & Chen (2008) introduced a similar theory in the context of stochastic dynamics of structures (Li & Chen 2009, ch. 7–8; Chen & Li 2009).

## 6. A numerical application to a tumour cell growth model

In order to verify the correctness of the aforementioned equation (5.6), we present here a numerical example. In particular, we consider the transient properties of the tumour cell growth model under immune response recently proposed by Zeng & Wang (2010). This model includes additive as well as multiplicative coloured noises (see also Fiasconaro *et al.* 2006) and it is described by the equations
6.1
where *x*(*t*;*ω*) denotes the population of tumor cells at time *t* while
6.2
In equation (6.2), *β* is the immune rate and *θ* is related to the rate of growth of cytotoxic cells. These parameters are typically set to *θ*=0.1, *β*=2.26. Also, the random process *f*(*t*;*ω*) represents the strength of the treatment (i.e. the dosage of the medicine in chemotherapy or the intensity of the ray in radiotherapy) while *η*(*t*;*ω*) is related to other factors, such as drugs and radiotherapy, that restrain the number of tumour cells. We shall assume that *f*(*t*;*ω*) and *η*(*t*;*ω*) are two independent Gaussians random processes with zero mean and correlation functions given by
6.3
where ℓ_{i} and *D*_{i} (*i*=1,2) denote, respectively, the correlation times and the correlation amplitudes^{6} of the processes *f*(*t*;*ω*) and *η*(*t*;*ω*). The factor 6 at the exponents has been introduced in order to make the correlation functions approximately zero when |*t*−*s*|≃ℓ_{i} (Venturi *et al.* 2008). Also, the initial condition *x*_{0}(*ω*) for the tumour density is set to be a standard Gaussian variable with mean 〈*x*_{0}(*ω*)〉=7.266. This mean value corresponds to the state of stable tumour (Zeng & Wang 2010) in the absence of random noises. We expand both processes *f*(*t*;*ω*) and *η*(*t*;*ω*) in a finite-dimensional Karhunen–Loève series
6.4
where {*ξ*_{1}(*ω*),…,*ξ*_{M1}(*ω*)} and {*ζ*_{1}(*ω*),…,*ζ*_{M2}(*ω*)} are two sets of zero-mean i.i.d Gaussian random variables, while *ϕ*_{k}(*t*) and *ψ*_{k}(*t*) are non-normalized eigenfunctions arising from the spectral decomposition of covariance kernels (6.3). The truncation of the series (6.4) is performed in order to retain 97 per cent of the total energy in the time interval [0,1]. The corresponding number of terms *M*_{1} and *M*_{2} is shown in table 1, as a function of the correlation time.

In figure 1, we plot several realizations of the response process *x*(*t*;*ω*) for different (randomly sampled) realizations of the forcing processes *f*(*t*;*ω*) and *η*(*t*;*ω*). The case with very small correlation times falls within the range of validity of the small correlation time approximation considered by Zeng & Wang (2010). The effective Fokker–Plank equation (see §4*b*) overcomes the curse of dimensionality, and it can be solved by standard numerical methods. However, for *large correlation times*, such approaches cannot be employed for obvious reasons.^{7} In these cases, we need to resort to other methods, e.g. methods based directly on random variables as discussed in §5. To this end, we consider the following joint response-excitation probability density associated with the system (6.1)–(6.4)
6.5
where the average is with respect to the joint probability density of the variables {*ξ*_{1}(*ω*),…,*ξ*_{M1}(*ω*)}, {*ζ*_{1}(*ω*),…,*ζ*_{M2}(*ω*)} and the initial condition *x*_{0}(*ω*). By differentiating equation (6.5) with respect to time and taking into account equation (6.1), we obtain
6.6
This is a linear transport equation in two variables, *t* and *a*, and *M*_{1}+*M*_{2} parameters. Under the assumption that the initial state of the system *x*_{0}(*ω*) is independent of {*ξ*_{1}(*ω*),…,*ξ*_{M1}(*ω*)} and {*ζ*_{1}(*ω*),…,*ζ*_{M2}(*ω*)}, the initial condition for the joint density (6.5) is explicitly given by
6.7
Once the solution to equation (6.6) is available, we obtain the response probability of the system by marginalizing (6.5) with respect to {*b*_{k}} and {*c*_{j}}, i.e.
6.8

The numerical solution to equations (6.1) and (6.6) is computed by using different approaches. Specifically, for the stochastic ODE problem (6.1), we have employed both Monte Carlo (5×10^{6} samples) and probabilistic collocation methods (Foo & Karniadakis 2008). In the latter case, depending on the number of random variables, we have used either Gauss–Hermite quadrature points or sparse grid (level 2) points for *ξ*_{k} and *ζ*_{j}. On the other hand, the partial differential equation (6.6) is first discretized in the *a* variable by using a Fourier–Galerkin spectral method of order 50, and then collocated at either Gauss–Hermite points or sparse grid points in the variables *b*_{k} and *c*_{j}. A second-order Runge–Kutta scheme is used to advance in time both equations (6.1) and (6.6).

In figure 2, we plot the time evolution of the response probability density of the system, i.e. the tumour density, for excitation processes with very large correlation time (ℓ_{i}=10) compared with the period of integration, which is [0,1]. We also show the effects of a variation in the correlation amplitude *D*_{2} characterizing the random process *η*(*t*;*ω*) restraining the number of tumour cells. The observed changes in the temporal evolution of the response probability density function are consistent with the results of Zeng & Wang (2010). The relevant statistics, i.e. the mean and the variance, of the tumour density are compared in figure 3 against similar results obtained by using different stochastic approaches. This comparison clearly shows that the transport equation (6.6) for the joint probability density function is indeed correct and allows for accurate predictions. The discrepancy observed in the lower right variance plot between the Monte Carlo solution (continuous line) and the sparse grid (level 2) solution of equation (6.6) (two variables and 11 parameters) is due to an insufficient integration accuracy when evaluating the statistical moments from the response probability density function.

## 7. Summary

By using functional integral methods, we have determined an equation describing the evolution of the joint response-excitation probability density function of a first-order nonlinear stochastic dynamical system driven by coloured random noise with arbitrary correlation time. This equation can be represented in terms of a superimposition of two differential constraints, i.e. two partial differential equations involving unusual limit partial derivatives. The first one of these constraints was determined by Sapsis & Athanassoulis (2008). We have addressed the question of computability of the joint response-excitation probability density function as a solution to a boundary value problem involving only one differential constraint. By means of a simple analytical example, we have shown that such a problem is undetermined, in the sense that it admits an infinite number of solutions. This result provides a definitive answer to the question first raised by Sapsis & Athanassoulis (2008, p. 295), regarding the solvability theory of the system of equations (2.8)–(2.10). In order to overcome this issue, we have included an additional differential constraint, i.e. the complementary constraint (2.15), which yields a complete evolution equation for the joint response-excitation density. This equation, however, involves an average requiring a closure approximation just like in the classical response approach (Moss & McClintock 1995). Such approximation can be constructed based on small correlation time expansions.

We have also studied nonlinear differential equations involving random parameters. This class of problems can be reformulated in an exact way, i.e. without any closure, in terms of a possibly high-dimensional linear transport equation for the joint response-excitation probability density function of the system. Such equation can be solved numerically by exploiting recent advances in numerical methods for high-dimensional systems such as proper generalized decomposition (Chinesta *et al.* 2010; Nouy 2010*a*), sparse grid collocation (Foo & Karniadakis 2008, 2010) or functional ANOVA techniques (Rabitz *et al.* 1999; Cao *et al.* 2009; Yang *et al.* submitted). In order to investigate this possibility, we have applied one of these methods, i.e. sparse grid collocation, to the evolution equation arising from the tumour cell growth model recently proposed by Zeng & Wang (2010). This allowed us to simulate the dynamics of the tumour density for Gaussian forcing processes with low to moderate correlation times, i.e. in a range where the small correlation time expansion considered by Zeng & Wang (2010) does not apply.

## Acknowledgements

We would like to thank Prof. B. L. Rozovsky of Brown University for stimulating discussions and very useful suggestions. This work was supported by OSD-MURI (grant no. FA9550-09-1-0613), by DOE (grant no. DE-FG02-07ER25818) and NSF (grant no. DMS-0915077).

## Appendix A. Differential constraints for the joint response-excitation probability density function of a nonlinear pendulum

In this appendix, we show that the functional integral approach leading to the differential constraint (2.8) can be extended, with minor modifications, to systems that can be more general than those described by equation (2.1). To this end, let us consider the dynamics of a nonlinear pendulum subject to an external random driving torque. A deterministic version of this problem has been extensively studied in the past as a prototype problem to understand chaos and the routes to chaos (D'Humieres *et al.* 1982; Blackburn *et al.* 1987; Gitterman 2005). The equation of motion is
A1
where *I* is the moment of inertia, *λ* the damping constant, the restoring torque, and *f*(*t*;*ω*) the random (coloured in time) driving torque. Equation (A1) can be written as a first-order system as follows:
A2
and
A3
The joint response-excitation probability density function for this system is
A4
Differentiation with respect to *t* yields
A5
Substituting (A1) into (A5) and taking the limit gives the following differential constraint for the joint probability density
A6
A partial integration of with respect to *ab* yields the marginal compatibility condition with the probability density function of the random torque acting on the oscillator
A7

## Appendix B. Small correlation time expansion for Gaussian forcing processes

In this appendix, we review a classical approach (Fox 1986) to determine a closure approximation to the averages appearing in equations (4.11) and (4.16) based on a small correlation time expansion.

(a) *Exponential correlation function*

Let us consider a Gaussian random forcing term with exponential correlation function
B1
where ℓ characterizes the *correlation time* of the process while *D* denotes the *correlation amplitude*. The parameter *D* ultimately controls the amplitude of the process *f*(*t*;*ω*). Let us calculate explicitly the last term appearing in equation (4.11). To this end, we first consider
B2
At this point, we expand the function
B3
in a power series around *s*=0. This is a key step for the approximation of the functional integral. Indeed, the variable *s* now represents the correlation time of the response process and we expect that for small ℓ we obtain small *s*. This is consistent with the white noise limit, where the response process is Markovian if the forcing process is white in time. Therefore, let us consider an expansion about the Markovian case
B4
For forcing processes with very small correlation time, we are allowed to retain only the linear term in the power series expansion (B4). This yields the following approximation of the term (B2)
and therefore
B5

(b) *Gaussian correlation function*

The approximation technique outlined in the previous subsection can be also applied to zero-mean Gaussian forcing processes having Gaussian correlation function B6 In this case, we find (for ) B7 where erf(⋅) denotes the standard error function.

## Footnotes

↵1 We shall assume that the process

*f*(*t*;*ω*) is real-valued and at least continuous with probability measure and sample space , which can be taken to be a quite general separable Banach space. For example, , , for . Standard existence and uniqueness theory (Sobczyk 2001) then ensures that there is a stochastic process*x*(*t*;*ω*) with sample space (i.e. at least differentiable), a probability measure and a joint probability space such that the joint process {*x*(*t*;*ω*),*f*(*t*;*ω*)} verifies the stochastic ordinary differential equation (2.1).↵2 The integral is formally written from to , although the probability density may be

*compactly supported*.↵3 In the numerical scheme employed by Sapsis & Athanassoulis (2008), this multiplicity was overcome by using a representation of the solution to equation (2.8) in terms of a Gaussian kernel. This introduces implicitly a symmetry in the covariance structure of the system, which ultimately results in a closure model.

↵4 In order to see this, let us recall that the cumulants

*C*_{n}(*t*_{1},…,*t*_{n}) of a Gaussian process are exactly zero for*n*≥3 and therefore equation (4.4), in this particular case, reduces to equation (4.5).↵5 We have assumed that the response process and the excitation process are independent at the initial time.

↵6 The correlation amplitude ultimately controls the amplitude of the process, namely, when

*D*_{1}is increased, then the amplitude of*f*(*t*;*ω*) increases.↵7 For large correlation times, the series representation (B4) cannot be truncated to first order in

*s*. This yields many additional terms in the closure approximation (B5), which ultimately render the closure itself not practical.

- Received March 23, 2011.
- Accepted October 17, 2011.

- This journal is © 2011 The Royal Society