Abstract
In this paper, we consider the problem of recovering the drift function of a Brownian motion from its distribution of first passage times, given a fixed starting position. Our approach uses the backward Kolmogorov equation for the probability density function (pdf) of first passage times. By taking Laplace transforms, we reduce the problem to calculating the coefficient function in a second-order ordinary differential equation (ODE). The inverse problem effectively amounts to finding the convection coefficient of the ODE, given the transformed pdf for positive values of the Laplace variable. Our first contribution is to find series solutions to the forward problem and show that the associated operator for the linearized inverse problem is compact. Our second contribution is numerical: for low noise levels, we reconstruct simple drift functions by applying Tikhonov regularization and performing a Newton iteration (Levenberg–Marquardt method). For larger noise, our solution displays large oscillations about the true drift.
Introduction
In this paper, we shall be concerned with finding the drift function of a Brownian process with constant diffusion coefficient, given the first passage time distribution (FPTD). Consider a particle that undergoes a random walk, , parametrised by a drift function and a constant diffusivity on the interval . Let the starting position be where . The particle cannot cross where there is a ‘reflecting’ boundary condition. When the value of first reaches 1, the time is recorded. Upon repeating the process many times, always starting at the same position , one obtains a distribution of first passage times. From this data, can one infer the drift function ?
We can get some intuition for this problem by considering the governing stochastic differential equation (SDE) explicitly,(1) (1) where , is the Wiener measure and is a constant diffusivity. In the case where , we have the deterministic equation . If the drift is always positive, the first passage time or exit time is given by(2) (2) and it is clear in this limit that from a single exit time, no pointwise information about can be gained: there are many functions that satisfy (Equation2(2) (2) ) for a given . If , we have a distribution of first passage times instead of a single time; now there is at least a possibility that some more information about can be extracted.
The study of these types of exit time problems is motivated by experiments that involve the dissociation of receptor–ligand complexes.[Citation1–Citation4] In these experiments, the force on a complex is increased until rupture occurs at a critical force. The experiment is repeated many times to obtain a distribution of rupture forces; thermal fluctuations ensure the rupture does not occur at a single fixed force. The objective of these experiments is to infer some quantitative features of the bond such as its dissociation constant. In the theoretical modelling of such systems, the rupture of the complex is assumed to proceed along a one-dimensional bond coordinate. Specifically, the problem can be treated as escape from a potential well.[Citation5–Citation8] Time dependent potentials (corresponding to the explicit application of a force ramp or a time-dependent drift function) can also be included.[Citation9] However, in this paper, we restrict ourselves to time independent potentials wells/drift functions.
One framework for approaching these stochastic inverse problems is to use the (deterministic) backward partial differential equation for the FPTD.[Citation10] After taking Laplace transforms, the problem amounts to calculating the coefficient function in a second-order differential equation. Mathematically, these types of problems have a long history, dating back to Borg.[Citation11] The aim of inverse Sturm–Liouville problems is to reconstruct the coefficient in the differential operator from its eigenvalues. Although spectra can be generated from the exit time distribution in certain cases, in our problem, we do not have direct access to eigenvalues. Nevertheless, a large body of work exists that discusses both analytic [Citation12, Citation13] and numerical methods [Citation14, Citation15] for coefficient reconstruction from its eigenvalues.
There are also several variations of the stochastic inverse problem where the quantity to be reconstructed and/or given data are different. We mention here a few recent theoretical and experimental studies. The assumption of constant in Equation (Equation1(1) (1) ) can be relaxed and could be inferred along with – this case was studied in [Citation16]. One could also infer the drift and diffusivity from knowledge of particle positions. Förster resonance techniques applied to ensembles of folding proteins [Citation17] give time snapshots of the probability distribution of from which it may be possible to infer . In [Citation18], functional forms of 2D drifts and a constant diffusivity were inferred from measurements of particle positions using Bayesian techniques.
In this paper, we we propose a Levenberg–Marquardt method to recover the drift function from the FPTD. The method is based on repeatedly solving a first-kind linear integral equation using regularization techniques. Unlike previous work, [Citation19] we do not a priori assume that the drift function can be decomposed into a pre-defined set of basis functions. At each stage, the regularization parameter is calculated using a method similar to Morozov’s discrepancy principle.
In Section 2, we derive the backward Kolmogorov equation that governs the exit time distribution in terms of a spatially dependent drift function. This equation defines a map from the drift to the exit time distribution. We also state the forward and inverse problems. In Section 3, we discuss our algorithm generally and prove compactness of the linearized map through five lemmas. In Section 4, we discuss more numerical aspects of our algorithm and the process of generating artificial data to test our method. Artificial data is required because in practice generating FPTDs by solving (Equation1(1) (1) ) was too time consuming; accurate approximation of the drift relied on having data with extremely low levels of noise. In Section 5, we present our numerical results. Finally, we summarize our findings with a conclusion in Section 6.
Statement of forward and inverse problems
Consider a Brownian particle that starts at a position and undergoes a Brownian motion parameterized by a spatially dependent drift function and unit diffusivity. We assume there is a hard wall at and if the particle exits at , the time of exit is recorded. When this process is repeated many times with the same starting position , we obtain a FPTD. In this section, we derive an equation for the FPTD in terms of the drift and then state the forward and inverse problems.
The probability that the particle is in the interval at time , given that it was at position at time is given by Kolmogorov’s forward equation [Citation10]:which is valid for . This equation is subject to the no flux and absorbing boundary conditions , and initial condition(3) (3) Kolmogorov’s backward equation can also be written for the same process (Equation1(1) (1) ), but in terms of initial position and time . The backward equation is formulated in terms of the adjoint operator :(4) (4) which is valid for . This equation is subject to the corresponding adjoint boundary conditions , and identical initial condition (Equation3(3) (3) ).
For time homogeneous processes, the probability distribution function only depends on the difference between final and initial times. Therefore, . Furthermore, we see that (Equation4(4) (4) ) is valid for and the time derivative can be rewritten asFinally, making the transformation (so that now ), we have(5) (5) The survival probability is the probability that the particle is ‘alive’ at time (the particle ‘dies’ when it exits) given that it started at position . The survival probability comes from integrating over all possible current positions : with since the particle is always alive at ; therefore the survival probability also obeys Equation (Equation4(4) (4) ). In fact, the Laplace transformed survival probability satisfies(6) (6) The exit time distribution, given that the particle started at position , can be calculated in terms of the survival probability . By definition, the probability of first exit in the time interval is given by . In this situation, the particle has survived until time but dies at time . Therefore, and . Their Laplace transforms are related through(7) (7) Upon multiplying both sides of (Equation6(6) (6) ) by and using (Equation9(9) (9) ), we have(8) (8) with corresponding boundary conditions derived from (Equation7(7) (7) ) and (Equation8(8) (8) )(9) (9) In the forward problem, we are given and the problem is to solve Equation (Equation10(10) (10) ) to find for regarding as a fixed parameter. Since the eigenvalues of (Equation10(10) (10) ), (Equation11(11) (11) ) and the homogeneous form of (Equation12(12) (12) ) are strictly negative, the solution to (Equation10(10) (10) )–(Equation12(12) (12) ) exists uniquely for any .
For the inverse problem, we assume that the Brownian process discussed above has been repeated many times from a starting point , yielding data for a fixed, known, and . Upon Laplace transforming, we obtain for . The inverse problem associated with (Equation10(10) (10) )–(Equation12(12) (12) ) is to find given . When , we have the trivial solution ; this simply means that , as expected since is a probability density function.
In the following sections, we will only refer to the Laplace-transformed exit time distribution . Therefore, we drop the tilde accent: symbols such as , and all refer to Laplace-transformed quantities.
Properties of the linearized map
Our numerical method is a Levenberg–Marquardt method, based on regularizing and then solving a first-kind integral equation at every iteration. We first describe the method in an abstract setting. Let be the forward map of (Equation10(10) (10) )–(Equation12(12) (12) ) so that . Most of this section is devoted to proving basic properties of the linearized map.
To solve the inverse problem for given data , we need to find a such thatThis equation can be linearized about a guess , so that small deviations satisfy(10) (10) where is the Frechet derivative of at (and hence changes with every iteration ). Below in Lemma 5, we show that is a compact operator and so has an unbounded inverse. Equations involving compact operators often give rise to ill-posed problems: therefore Equation (Equation13(13) (13) ) may have multiple solutions or no solution at all. One way of remedying this problem is to find a that minimizes a Tikhonov functional:(11) (11) where is a small regularization parameter. The second integral represents a penalty term which ensures that cannot be too large. For noiseless data , when in (Equation14(14) (14) ), it can be shown [Citation20] that converges to the least squares solution to (Equation13(13) (13) ). Furthermore, the minimizer of (Equation14(14) (14) ) is unique in satisfying [Citation21](12) (12) where is the identity operator and is the adjoint of . Deciding on a ‘good’ choice of in (Equation15(15) (15) ) at each step of the iteration is an open problem [Citation22]; in our implementation, given the st iterate , note that is a function of . Its value is chosen to minimize where is a (non-uniform) discretization of . Further details are discussed in Section 4. The numerical solution of (Equation15(15) (15) ) forms the backbone of our algorithm. Solving for the Newton iterates in (Equation13(13) (13) ) using (Equation15(15) (15) ), along with a method for finding , constitutes the Levenberg–Marquardt algorithm. Some convergence properties of this algorithm have been proved by Hanke in [Citation23].
For concreteness, we now find the explicit form of . Linearizing (Equation10(10) (10) ) about a known solution , small changes are related to small changes through(13) (13) Since both and must obey the boundary conditions (Equation11(11) (11) ) and (Equation12(12) (12) ), we have and . Given for fixed and , we compute by solving the first-kind integral equation(14) (14) where is the Greens function for Equation (Equation16(16) (16) ).
In summary, given a Laplace-transformed FPTD and a known , we first make a guess for the drift, . Using , we solve the forward problem Equations (Equation10(10) (10) )–(Equation12(12) (12) ) to find and compute for . Regularizing and then solving (Equation17(17) (17) ) then yields a which we use to update our guess. To implement our algorithm, we need the form of the Greens function in Equation (Equation17(17) (17) ). This is provided through the following lemma.
Lemma 1
Let . The Greens function that satisfies(15) (15) has the form(16) (16) where satisfies (Equation10(10) (10) )–(Equation12(12) (12) ) and the second linearly independent solution satisfies the same ordinary differential equation (Equation10(10) (10) ) but with modified boundary conditions and . Moreover, exists (is bounded) for .
Proof
By writing (Equation10(10) (10) ) in self adjoint form, we deduce that the Greens function must satisfywhich is identical to Equation (Equation18(18) (18) ). Let be a solution to (Equation10(10) (10) ) that satisfies , and is linearly independent from . Then we can write asand satisfies the boundary conditions (Equation19(19) (19) ) and (Equation20(20) (20) ). In addition, at , must be continuous at and its derivative must jump by . These conditions yield(17) (17) where . Finally, by Abel’s identity,(18) (18) for some ; evaluating (Equation25(25) (25) ) at gives . Hence (Equation23(23) (23) ) and (Equation24(24) (24) ) imply and .
Because , in Equation (Equation22(22) (22) ) always exist for and . Furthermore, is always non-zero: if , the boundary condition (Equation11(11) (11) ) along with (Equation10(10) (10) ) implies so that the right boundary condition (Equation12(12) (12) ) cannot be satisfied. Therefore, is bounded for and . (We prove that is also finite at as part of Lemma 4.)
In the case where , we have an explicit form for the Green’s function, whose properties are used extensively in Lemmas 3, 4 and 5.
Lemma 2
The Green’s function satisfyingfor has the following properties:
.
For all , .
For all , .
Proof
(a) Solving explicitly, we have(19) (19) (20) (20) and follows immediately.
(b) and (c): These results follow immediately from Equations (Equation26(26) (26) ) and (Equation27(27) (27) ) and using , and for .
Lemma 3
Let for all . Then a series solution to (Equation10(10) (10) )–(Equation12(12) (12) ) exists in the form(21) (21) where(22) (22) for with the convention for any sequence . The Green’s function is defined in (Equation26(26) (26) ), prime denote differentiation with respect to the function’s first argument and .
Proof
Using (Equation26(26) (26) ), we can rewrite (Equation10(10) (10) ) in terms of an integro-differential equation(23) (23) Substituting (Equation28(28) (28) ) into (Equation31(31) (31) ), we find that and for . These relations imply (Equation29(29) (29) )–(Equation30(30) (30) ). To finish the proof, we now show that the series (Equation28(28) (28) ) is convergent. From Lemma 2(b) and (c), we haveTherefore,and the partial sum is bounded from above and increasing. It follows that converges and therefore so does .
In light of Lemma 3, it is now straightforward to show that the inverse problem (Equation10(10) (10) )–(Equation12(12) (12) ) does not admit a unique solution . This is proved in the next two corollaries. Compactness of the linearized map is established in Lemmas 4 and 5.
Definition 1
For any function , fix and define as(24) (24) Note that is not necessarily continuous even if is. However is always piecewise continuous in if is continuous: see Figure .
Corollary 1
Assume that there exists a constant such that . Fix the starting position . Then there is a such that is invariant if the drift in Equation (Equation10(10) (10) ) is replaced by . Therefore, when , the two drift coefficients and give identical for any and therefore identical FPTDs for all .
Proof
Let . We show that (Equation29(29) (29) )–(Equation30(30) (30) ), and therefore, (Equation28(28) (28) ) is invariant when is replaced by . We partition the range of integration into two setsso that for ,where(25) (25) and . Clearly, is invariant when is replaced with since when . In where for , we substitute (so for ) and immediately obtain . The terms in the integrand then transform as follows:(26) (26) using (Equation26(26) (26) ). This equation holds if and .
The products of derivatives of satisfy
Finally, products of satisfy
For a general , the transformation will yield a that is discontinuous (see Figure ) and therefore inadmissible – recall we are looking for continuous to the inverse problem (Equation10(10) (10) )–(Equation12(12) (12) ). However, there are also infinitely many drift functions where and are both continuous as illustrated in Figure and in the next corollary. In the extreme case where , we can take : and are both continuous and generate identical FPTDs.
Corollary 1 complements the results in [Citation16]. Theorem 1 in [Citation16] implies that if the diffusivity in the Brownian motion is a known constant and the drift is known on where (using the notation in [Citation16]), then knowledge of uniquely determines on , where the measurement point and is a countable set of measure zero in (0,1). This does not contradict corollary 1. For , the authors in [Citation16] essentially proved that this a priori knowledge of on resolves the vs. ambiguity. On the other hand, corollary 1 proves that when , the a priori knowledge is insufficient to determine uniquely. For example, we can take and and use corollary 1 to construct different and that are identical for . One could also consider the case where the measurement point , is known to the left of the measurement point but unknown to the right, corresponding to Theorem 2 in [Citation16]. Then the theorem states that a sufficient condition for to be uniquely determined in is that . Corollary 1 shows that this condition is also necessary.
Corollary 2
Fix the starting position and fix . Let be a continuous function such that and when . Let be the Laplace-transformed FPTD for . Then if is perturbed by , and have identical FPTDs.
Proof
and generate the same and .
Before we show compactness of the linearized operator , we need an intermediate result. We prove
Lemma 4
Let for and fix . Thenwhere is defined in Lemma 1 and satisfies (Equation10(10) (10) )–(Equation12(12) (12) ).
Proof
We separate the proof into three steps. Below, we let and make frequent use of the inequalities , and for .
First, we haveFor ,with the convention for any sequence . Therefore(27) (27)
Second, the Greens function obeys (Equation18(18) (18) ) which we write as(28) (28) We can write the solution to (Equation38(38) (38) ) in terms of the Greens function for , :(29) (29) We look for a series solution to this integro-differential equation in the form . We follow a similar procedure to the Proof of Lemma 3 and substitute this series into (Equation39(39) (39) ) to findfor , where and prime denotes differentiation with respect to the function’s first argument. Using Lemma 2(b) and (c) and (see the definition in Equation (Equation21(21) (21) )),for . Therefore(30) (30)
Finally, using (Equation37(37) (37) ) and (Equation40(40) (40) ),(31) (31)
We finally arrive at a compactness result for (see Equation (Equation13(13) (13) )). This property of means that is unbounded and justifies the use of Tikhonov regularization via Equation (Equation14(14) (14) ).
Lemma 5
Let for all . Fix and . Define by(32) (32) where is defined by (Equation22(22) (22) ) and satisfies (Equation10(10) (10) )–(Equation12(12) (12) ). Then maps bounded subsets of to relatively compact subsets in . Hence is a compact operator from to .
Proof
First, we show that maps functions from to functions in . Using Schwarz’s inequality, we see that(33) (33) where is defined in Equation (Equation42(42) (42) ). This function is square integrable since(34) (34) Now, we show that maps bounded subsets to relatively compact subsets. Let be a bounded subset of and let . Then there is a constant such that(35) (35) By showing that is bounded and equicontinuous, we prove that it is relatively compact through the Frechet–Kolmogorov theorem [Citation24], an extension of the Arzela–Ascoli theorem for for spaces. (Recall that is equicontinuous in if , such that .)
is bounded: From (Equation44(44) (44) ) and (Equation45(45) (45) ), . Hence is bounded in .
is equicontinuous in : Note that we have(36) (36) which is proved in the appendix; see (58). Classical solutions to (Equation10(10) (10) )–(Equation12(12) (12) ) for are differentiable so that . The Greens function is defined in terms of classical solutions and in (Equation22(22) (22) ) so that we also have . Therefore . Since this function is bounded at by (Equation41(41) (41) ), it is also uniformly continuous in . Uniform continuity of implies that , (depending only on and possibly ) such that , and ,(37) (37) Furthermore, we have(38) (38) for . Inequalities (Equation41(41) (41) ) and (Equation47(47) (47) )–(49) imply(39) (39) uniformly in . Now choose small enough so that . It follows that for each , , (depending on and , but not on and ) such that and ,Therefore, , and , independent of and such that impliesNow fix and define . Choose so that , with defined through (Equation46(46) (46) ). Then for all ,where we define whenever the argument is negative. Hence is equicontinuous in .
Algorithm for drift reconstruction
Lemma 5 shows that the operator is compact and so locally, the inverse problem (Equation10(10) (10) )–(Equation12(12) (12) ) is ill-posed: there may be no solution for , more than one solution or the solution may not vary continuously with . In fact, we showed a stronger result in Corollary 1: that globally the inverse problem (Equation10(10) (10) )–(Equation12(12) (12) ) is not unique. Here, we implement a numerical method based on regularizing and solving the first-kind integral equation (Equation17(17) (17) ). The regularization of the integral equation alleviates the issue of local non-uniqueness. Global non-uniqueness can be partially remedied by only using Laplace-transformed FPTDs that correspond to symmetric drifts () as data to the inverse problem (Equation10(10) (10) )–(Equation12(12) (12) ).
We now describe the main steps of the Levenberg–Marquardt algorithm. We represent the current guess for our solution and the update on a chebyshev grid , so thatIn addition, for a fixed integer , letThe evaluation nodes are always finite.
Given a fixed , noisy exit time data (generated using the method described in Section 4.1), a trial drift function a tolerance level, and the maximum number of iterations allowed, , the algorithm to reconstruct from is as follows:
Let . While and :
Solve Equation (Equation10(10) (10) ) with and , to compute the th iterate . We use a pseudospectral method based on the cheb.m routine in [Citation25].
Compute .
Find that minimizes(40) (40) where is the numerical solution to (Equation10(10) (10) )–(Equation12(12) (12) ) with drift functionwith satisfying(41) (41) and(42) (42) In (52), and is the identity. In Equation (53) are weights for Clenshaw–Curtis quadrature on and are the abscissae. These weights and abscissae are related to the ones on [1,1], , through , . Matlab can quickly generate using routines such clencurt.m.[Citation25] The sum (51) was minimized using Matlab’s fminbnd.m with a lower bound for of and an upper bound of .
Set .
Let .
Minimizing (51) amounts to selecting using a variation of Morozov’s discrepancy principle. We also tried to determine by minimizing instead of (51), with appropriate rescaling of the infinite integration range, but this method generally resulted in a scheme with worse convergence properties.
Generation of noisy data
The Laplace transform of perfect first passage time data corresponding to a target drift is generated by solving Equation (Equation10(10) (10) ) for . Using , we generate noisy data through(43) (43) where tilde denotes Laplace transform and is the Laplace transform of a ‘noisy’ function. We use Gaussian noise for discretized as which are generated through(44) (44) for . where , and . Hence, are normally distributed with mean 0 and standard deviation ; controls the frequency of the noise and is the cut-off beyond which . The statistics of the noise are completely specified through the three parameters .
Given a uniform-in-time discretization its Laplace transform at is numerically implemented through(45) (45) where . In Equation (56), we are approximating the Laplace integral using a one-sided rectangle rule.
One can think of the as being the difference between the exact FPTD at times and a histogram of simulated first passage times using bins; however, Equation (54) represents the noise of the exit time distribution only qualitatively. We test our algorithm against noisy data to prevent an ‘inverse crime’; a natural extension is to generate first passage times by numerically integrating (Equation1(1) (1) ) in time. In particular, since for , we have assumed that the noise becomes insignificant providing the first passage times are sufficiently large. In reality, noise will persist for large but it will be small in magnitude: for a finite number of realizations, there will be a maximal exit time; beyond this time, the numerical FPTD is exactly zero. Since is a probability density and must be integrable on , as . Their difference (which represents the noise) as .
Results and discussion
Some recovered drift functions are shown in Figure . In these plots, we test our method with noisy data and vary the noise parameters in Equation (55) to test for convergence. In (a) and (b), we vary the standard deviation (amplitude) of the noise, . We see that the noise has to be very small for the method to converge. When is small enough, our method can capture the general shape of some simple drift functions. We point out that since the added noise is random, two data-sets with the same can give different results; for example, although reasonable results are obtained with in (a), the method can also diverge for a different noise realization with the same . In both (a) and (b), when the noise level becomes large enough ( for these cases), the method produces oscillations about the true drift. These oscillations start small and their growth saturates after many iterations: the oscillatory solutions in Figure are all ‘steady state’. To illustrate the sensitivity of the method to noise, we point out that in (a) and (b) when , and , the reconstructed drifts contain large errors and yet for these noise parameters, the (Laplace transformed) noisy data only differ from the (Laplace transformed) perfect data by .
In (c), we vary the timespan of the added noise. We see that larger values of render the method unstable. This suggests that errors in the largest first passage times are more detrimental to drift reconstruction than smaller ones. The implication is that to find in practice, one has to perform enough realizations of the Brownian motion so that the tails of the FPTD are captured accurately. Finally in (d), we vary the noise frequency which in practice could be related to the number of bins in the histogram. Because low frequency noise gives poorer results than high frequency noise, increasing the number of bins (while keeping the noise amplitude constant) in the numerical FPTD could give more accurate reconstructions. Recall that representation of as a polynomial is inherent in our numerical method. Because the target drift in (d) is not differentiable, it is not surprising that large oscillations are always present, even for noiseless data.
In Figure we plot the two-norm of the difference between the exact solution and the iterates , , as a function of iteration number, . For small values of , the error in (a) after about iterations is still on the order of and does not seem to decrease. Although the convergence properties of our method still need to be further explored, existing studies on regularized Gauss–Newton methods in the context of inverse scattering [Citation26] suggest that convergence of such methods can be logarithmic in . As increases, (a) shows that the method settles to a solution that is further away from the target. Also, it is more likely that large oscillations develop in the numerical solution after a promising start, as shown in the case. This phenomenon of semi-convergence [Citation22] with noisy data is common in inverse problems and can be partially remedied through stopping rules.[Citation22, Citation27] In (b), we plot the values of the regularization parameters as a function of the iteration number (there is no regularization for where the iterate is just the initial guess). At the th iteration, the shown in (b) are the calculated regularization parameters used to find in (a). It appears that the best approximation to is obtained when , whose value is determined by the minimization problem (51), rapidly decreases with .
In Equation (Equation17(17) (17) ), we see that since when and , rows and of in (53) consist only of zeros, and columns and of consists only of zeros. Therefore, (52) implies that(46) (46) and boundary values of the iterates do not change. In Figure the exact solution has and our method fails to converge to the exact solution even for perfectdata. Instead, because the initial guess is , the end points of the solution are pinned at zero and the rest of the solution oscillates about the target . Like all Newton methods, the choice of initial guess is vital to the success of our algorithm and this is illustrated in (b). For initial guesses that have , the method is unable to reconstruct even basic features of . However, when , the method quickly converges to the target drift.
We did attempt to reconstruct drift functions from noisy FPTDs generated by numerically integrating the SDE (Equation1(1) (1) ) using a basic Euler–Maruyama method. However, there are(at least) two sources of error associated with sampling in this way. Thefirst is related to the finite discretization of the SDE; this method also over predictsthe first passage times.[Citation28] The second is due to the finite sample size. When we used and 4,000 first passage times (which took about 1.5 hours to generate on a quad core Intel Xeon 2 GHz workstation), the error in was about 0.005 which waslarge enough to prevent convergence of the Levenberg–Marquardt method: ourmethod is extremely sensitive to small amounts of noise in the Laplace-transformeddata.
In Figure , data with usually produced convergent results; these noise parameters for the artificial data correspond to in Equation (54). Unfortunately, we found that it was extremely time consuming to generate a FPTD through the SDE whose Laplace transform is accurate to . Since the error in the Laplace-transformed data scales as the inverse square-root of the sample size, while the time taken to generate the data scales linearly with the sample size, we estimate that it would take about samples (with ) to obtain FPTD data that is accurate to . Therefore, generating samples would take about 40 years on our workstation. We leave as future work the implementation of more advanced methods [Citation28, Citation29] that can simulate diffusive FPTDs more efficiently.
Conclusions
In this paper, we studied the reconstruction of a drift function in a Brownian motion from the Laplace transform of its FPTD. Our main contributions are a proof that the linearized inverse problem is ill-posed and a numerical method based on the Levenberg–Marquardt method. This method allows drift reconstruction to be implemented from artificial data, generated by adding Gaussian noise to a solution of the forward problem.
Our numerical method is able to reconstruct simple drift functions for low noise levels. Furthermore, as with many Newton-type methods, convergence is contingent on a suitable initial guess. In particular, unless the initial guess already takes the correct values at the endpoints, the iterates do not converge.
We see our future work as consisting of three main parts. First, the inverse problem should be solved using exit time data generated by the SDE (Equation1(1) (1) ). From our results so far, we anticipate that the Brownian dynamics must be simulated quickly and accurately in order to generate exit time distributions with sufficiently low levels of noise. Alternatively, other efficient methods to sample from the FPTD must be developed.
Second, it remains to establish whether the assumption of is sufficient to guarantee uniqueness of the inverse problem (Equation10(10) (10) )–(Equation12(12) (12) ). When is not symmetric, we showed in corollary 1 that is not unique. In particular, if the starting position, and generate exactly the same distribution of first passagetimes.
Finally, our numerical method should be extended so that it can recover drift functions that have non-zero boundary values. An important first step would be to establish the behaviour of in (Equation43(43) (43) ) as .
Acknowledgments
The author thanks Fioralba Cakoni and Tom Chou for helpful discussions. This work was supported by the University of Delaware Research Foundation (UDRF).
References
- Fuhrmann A, Anselmetti D, Ros R, Getfert S, Reimann P. Refined procedure of evaluating experimental single-molecule force spectroscopy data. Phys. Rev. E. 2008;77:031912–1-31912-10.
- GetfertS, EvstigneevM, ReimannP. Single-molecule force spectroscopy: practical limitations beyond Bell’s model. Physica A. 2009;388:1120–1132.
- Getfert S, Reimann P. Suppression of thermally activated escape by heating. Phys. Rev. E. 2009;80:030101–1-030101-4.
- HeymannB, GrubmüllerH. Dynamic force spectroscopy of molecular adhesion bonds. Phys. Rev. Lett. 2000;84:6126–6129.
- BalseraM, StepaniantsS, IzrailevS, OonoY, SchultenK. Reconstructing potential energy functions from simulated force-induced unbinding processes. Biophys. J. 1997;73:1281–1287.
- DudkoOK. Single-molecule mechanics: new insights from the escape-over-a-barrier problem. Proc. Nat. Acad. Sci. U.S.A. 2009;106:8795–8796.
- DudkoOK, HummerG, SzaboA. Theory, analysis, and interpretation of single-molecule force spectroscopy experiments. Proc. Nat. Acad. Sci. U.S.A. 2008;105:15755–15760.
- HummerG, SzaboA. Kinetics from nonequilibrium single-molecule pulling experiments. Biophys. J. 2003;85:5–15.
- FreundLB. Characterizing the resistance generated by a molecular bond as it is forcibly separated. Proc. Nat. Acad. Sci. U.S.A. 2009;22:8818–8823.
- GardinerCW. Handbook of stochastic methods. Berlin: Springer; 1985.
- BorgG. Eine Umkehrung der Sturm-Liouvilleschen Eigenwertaufgabe. Acta Math. 1946;78:1–96.
- McLaughlinJR. Analytical methods for recovering coefficients in differential equations from spectral data. SIAM Rev. 1986;28:53–72.
- BrownBM, SamkoVS, KnowlesIW, MarlettaM. Inverse spectral problem for the Sturm-Liouville equation. Inverse Probl. 2003;19:235–252.
- RundellW, SacksPE. The reconstruction of Sturm-Liouville operators. Inverse Probl. 1992;8:457–482.
- Zhou S, Cui M. Determination of unknown coefficients in parabolic equations. J. Heat Transfer. 2009;131:111303–1-111303-6.
- BalG, ChouT. On the reconstruction of diffusions from first-exit time distributions. Inverse Probl. 2003;20:1053–1065.
- LipmanEA, SchulerB, BakajinO, EatonWA. Single-molecule measurement of protein folding kinetics. Science. 2003;301:1233–1235.
- MassonJ-B, CasanovaD, TürkcanS, VoisinneG, PopoffM, VergassolaM, AlexandrouA. Inferring maps of forces inside cell membrane microdomains. Phys. Rev. Lett. 2009;102:1–4.
- Fok P-W, Chou T. Reconstruction of potential energy profiles from multiple rupture time distributions. Proc. R. Soc. London, Ser. A. 2010;466:3479–3499.
- Groetsch CW. Integral equations of the first kind, inverse problems and regularization: a crash course. J. Phys.: Conf. Ser. 2007;73:45–54.
- KressR. Linear integral equations. Berlin: Springer-Verlag; 1989.
- EnglHW, HankeM, NeubauerA. Regularization of inverse problems. Dordrecht: Kluwer; 1996.
- HankeM. A regularizing Levenberg-Marquardt scheme with applications to inverse groundwater filtration problems. Inverse Probl. 1997;13:79.
- ChipotM. Elliptic equations: an introductory course. Basel: Birkhäuser Verlag; 2009.
- TrefethenLN. Spectral methods in Matlab. Philadelphia (PA): SIAM; 2000.
- HohageT. Logarithmic convergence rates of the iteratively regularized Gauss-Newton method for an inverse potential and an inverse scattering problem. Inverse Probl. 1997;13:1279.
- Colton D, Engl HW, Louis AK, McLaughlin JR, Rundell W, editors. Surveys on solution methods for inverse problems. New York: Springer; 2000.
- BuchmannFM. Simulation of stopped diffusions. J. Comput. Phys. 2005;202:446–462.
- GiraudoMT, SacerdoteL. An improved technique for the simulation of first passage times for diffusion processes. Comm. Statist. Simulation Comput. 1999;28:1135–1163.
Inequality for uniform continuity
In this appendix, we prove the inequality (Equation47(47) (47) ). Let where for ease of notation, we omit the and dependence in . Then for ,Therefore,(47) (47)