Jump Locations of Jump-Diffusion Processes with State-Dependent Rates

TL;DR

This paper introduces a framework for analyzing jump locations in jump-diffusion processes with state-dependent jump rates, revealing their relationship with the overall process distribution and demonstrating applications with non-monotonic behaviors.

Contribution

It develops an iterative map approach to compute jump location distributions in state-dependent jump-diffusions, linking jump locations to the stationary density of the process.

Findings

Derived an iterative method for jump location distributions.

Linked long-term jump location distribution to stationary density.

Showed non-monotonic dependence of jump behavior on diffusion strength.

Abstract

We propose a general framework for studying jump-diffusion systems driven by both Gaussian noise and a jump process with state-dependent intensity. Of particular natural interest are the jump locations: the system evaluated at the jump times. However, the state-dependence of the jump rate provides direct coupling between the diffusion and jump components, making disentangling the two to study individually difficult. We provide an iterative map formulation of the sequence of distributions of jump locations. Computation of these distributions allows for the extraction of the interjump time statistics. These quantities reveal a relationship between the long-time distribution of jump location and the stationary density of the full process. We provide a few examples to demonstrate the analytical and numerical tools stemming from the results proposed in the paper, including an application that shows a non-monotonic dependence on the strength of diffusion.

Equations

$$\tilde{p}(x,t)\text{d}x:=\mathsf{P}\left[X_t\in (x,x+\text{d}x)\right].$$

$${\partial }_t\tilde{p}(x,t)=\mathbb{L}\tilde{p}-\lambda (x)\tilde{p}+\mathbb{J}\lambda \tilde{p},$$

$$\text{d}{Y}_t=A(Y_t)\text{d}t+\sqrt{2D(Y_t)}\text{d}{W}_t,$$

$${\partial }_{t}p(y,t)=\mathbb{L}p:=-{\partial }_y\left\{A(y)p\right\}+{\partial }_{yy}\left\{D(y)p\right\}.$$

$$X_{{\tau }_{+}}=X_{{\tau }_{-}}+\Delta .$$

$$\mathbb{J}p(x,t):=p(x-\Delta ,t).$$

$$X_{{\tau }_{+}}=\eta .$$

$$\mathbb{J}p(x,t):=\delta (x-\eta ){\int }_{-\infty }^{\infty }p(x,t)\text{d}x,$$

$$\mathbb{J}p(x,t):={\int }_{-\infty }^{\infty }p(x-\Delta ,t)\mu (\Delta ,x)\text{d}\Delta .$$

$$X_{{\tau }_{+}}=\gamma X_{{\tau }_{-}}⟹\mathbb{J}p(x,t):=\gamma p(\gamma x,t)$$

$$\{\begin{array}{l}{\partial }_{t}p(x,t)=\mathbb{L}p-\lambda p\\ {\partial }_{t}q(x,t)=\lambda p.\end{array}$$

$$p_{\ell }(x)=q(x,\infty )={\int }_0^{\infty }\lambda (x)p(x,t)\text{d}t$$

$$p_{\tau }(t)={\int }_{-\infty }^{\infty }\lambda (x)p(x,t)\text{d}x.$$

$$s(t):=\mathbb{P}[\tau >t]={\int }_{-\infty }^{\infty }p(x,t)\text{d}x,$$

$$p_{\tau }(t)=-\frac{\text{d}s}{\text{d}t}=-{\int }_{-\infty }^{\infty }{\partial }_{t}p\text{d}x.$$

$$p_{\tau }(t)={\int }_{\infty }^{\infty }\lambda p\text{d}x={\int }_{\infty }^{\infty }{\partial }_{t}q\text{d}x.$$

$$p_{\ell }(x)={\int }_0^{\infty }{\partial }_{t}q\text{d}t=q(x,\infty )-q(x,0)=q(x,\infty )$$

$$p_i(x)\text{d}x:=\mathsf{P}[X_t\in (x,x+\text{d}x)|t=t_i].$$

$$\{\begin{array}{l}{\partial }_t{\hat{p}}_i(x,t)=\mathbb{L}{\hat{p}}_i-\lambda {\hat{p}}_i\\ {\hat{p}}_i(x,0)=\{\begin{array}{ll}\mathbb{J}{p}_{i-1}& i>1\\ p_0& i=1\end{array}\\ p_{i+1}(x)={\int }_0^{\infty }\lambda {\hat{p}}_i\text{d}t,\end{array}$$

$${\hat{p}}_i(x,t)\text{d}x:=\mathsf{P}[X_{t+t_{i-1}}\in (x,x+\text{d}x),t+t_{i-1}<t_i].$$

$$u_i:=p_i(x)/\lambda (x).$$

$$\mathbb{T}{u}_{i+1}:=\left[\lambda (x)-\mathbb{L}\right]u_{i+1}=\mathbb{J}\lambda u_i,$$

$$u_{i+1}={\mathbb{T}}^{-1}\mathbb{J}\lambda u_i.$$

$${\hat{p}}_i(x,0)={\int }_0^{\infty }\mathbb{L}{\hat{p}}_i\text{d}t-\lambda (x){\int }_0^{\infty }{\hat{p}}_i\text{d}t.$$

$$u_1(x)\to u_2(x)\to \cdots \to u_{\star }(x)\to u_{\star }(x)\to \cdots ,$$

$$0=\mathbb{L}{u}_{\star }-\lambda u_{\star }+\mathbb{J}\lambda u_{\star },$$

$$0=\mathbb{L}{\hat{p}}_s-\lambda {\hat{p}}_s+\mathbb{J}\lambda {\hat{p}}_s.$$

$${\int }_{-\infty }^{\infty }{\hat{p}}_s\text{d}x=1.$$

$${\int }_{-\infty }^{\infty }{u}_{\star }(x)\lambda (x)\text{d}x=1.$$

$$⟨{\tau }_i⟩={\int }_{-\infty }^{\infty }{u}_i\text{d}x={\int }_{-\infty }^{\infty }\frac{p_i(x)}{\lambda (x)}\text{d}x.$$

Full text

1 Introduction

Stochastic models driven by both Brownian noise (diffusion) and a discrete jump component, or so-called jump-diffusion models, have been used to describe a wide variety of phenomena. Perhaps most prominently, jump-diffusion has seen widespread use in mathematical finance [57, 36, 37, 24, 5] to describe frequent small transactions with occasional larger movements. Jump-diffusion has also served useful in mathematical biology, describing the integrate-and-fire nature of neuronal dynamics [28, 51, 16]. Other applications include biophysical descriptions of movement of chromosomes [52], interaction between soil moisture and rainfall events [12], and the occurrence of radio pulsar glitches [26]. Also relevant to the results described in this paper is a subset of jump-diffusion models that neglect diffusive noise and are driven by deterministic dynamics between jumps. This type of model appears in an equally eclectic variety of applications [54, 14, 18, 4].

Of particular interest in this work are jump-diffusion processes with state-dependent jump rates. This subset is seen in all the aforementioned applications, including finance [5, 29], ecology [18], biology [13], and astronomy [26]. This is often a natural supposition to make of a model, as the rate at which the jump event occurs may not remain constant on the timescale of interest, but instead depend on the state of the system. For instance, a neuron firing is dependent on the current voltage or a financial asset is often more likely to crash as it rises in price.

For many of these applications, the inclusion of a jump component means that it is a significant feature to the modeler. Thus, it is natural to ask: when (in some sense) do these jumps occur? While this often leads to studying the interjump times, a different possible interpretation (and the focus of this work) is to study the jump locations: the system evaluated at the jump times. This could be the voltage at which the neuron fires (the “threshold voltage”) or the price at which a stock crashes. Because the jump intensity is state-dependent, and the state itself is governed by a stochastic process, these locations are inherently random and difficult to disentangle from the full process. This raises another question: in what sense is knowledge (e.g. statistics) of the full process equivalent to knowledge of the jump locations? Does knowledge of the jump locations provide other insight to the behavior of the system? In this work, we seek to investigate these issues.

Previous studies have thoroughly investigated jump-diffusion models for constant jump rates [12, 30], as the two sources of noise can be independently superimposed. The interjump time distributions for processes with no diffusive noise have also been studied [13, 14] primarily using renewal theory or exploiting the deterministic nature between jumps. However, in general, state-dependent jump-diffusion is not deterministic between jumps, nor a renewal process, as sequential distributions are not independent, hence previous approaches do not apply. The theory of Cox processes [11] is also well-established and refers to a diffusion process driving a state-dependent Poisson process, but this does not complete the feedback loop of the jump process further modifying the state, which our framework allows. Diffusion with switching, which bears a resemblance to jump-diffusion, has also been studied [41] but is distinctly different in the role of discrete noise and not addressed in this work. Although efficient Monte Carlo methods exist for jump-diffusion and Lévy processes [56, 29], the jump component may be rare and therefore expensive to find an accurate distribution empirically. The solution to the full PDE describing this process can also be computed, but it is seemingly unknown how to disentangle the jump component alone from this description.

In this work, we present a general (accommodating a wide variety of models) framework for studying jump-diffusion systems and focus particularly on studying the sequence of jump locations. We present an iterative map formulation explicitly describing the distribution of the $i$th jump location. Fron this sequence of jump locations, statistics of the interjump times can be extracted even with diffusion, overcoming a limitation of previous works [14]. By taking the limit of the map, assuming the process reaches stationery, we find that the density of the full process differs from that of the jump locations if and only if the jump rate is state-dependent. An explicit relationship between these two distributions is established, meaning that knowledge of one immediately provides knowledge of the other.

A few examples are discussed, illustrating when this framework can be used to elucidate interesting behavior of a system. The first example, a stochastic model of neuronal integrate-and-fire, is used to demonstrate that the jump locations themselves (in this case, the firing voltages) may be interesting. For the second example, we propose a simple model of a molecular motor (such as kinesin [31]) taking force-dependent steps pulling a diffusive cargo. Because the stepping rate is state dependent, but steps (and cargo diffusion) also modify the state, the effective stepping rate (a proxy for the motor’s ability to produce transport) is non-trivial and can be studied using the proposed framework.

2 Formulation

2.1 Jump-diffusion

Let $X_{t}$ denote a state-dependent jump-diffusion process and $\tilde{p}(x,t)$ be the probability density (PDF) of the full process, defined formally as

[TABLE]

The evolution of the density is described by the forward Champman-Kolmorgov equation [27, 7]

[TABLE]

where the operator $\mathbb{L}$ characterizes the diffusion component, $\lambda(x)$ is the state-dependent intensity at which the jump process occurs, and the operator $\mathbb{J}$ describes the jump. A feature we seek to emphasize is the state dependence of the jump-rate $\lambda(x)$, resulting in direct coupling between the diffusion and jump components of the process.

The particular choices of $\mathbb{L},\mathbb{J},\lambda$ significantly change the behavior and characteristics of the process, but the description (2) is flexible enough to accommodate a wide variety of models, upon which we now elaborate.

2.2 Diffusion component

Consider the classical one-dimensional Itô SDE driven by Brownian motion described by

[TABLE]

where $W_{t}$ is a Brownian motion. The Fokker-Planck (forward Chapman-Kolmogorov) equation describing the evolution of the probability density function $p(y,t)$ corresponding to the path-wise description (3) is

[TABLE]

We use the Fokker-Planck operator $\mathbb{L}$ to characterize the diffusion process for the remainder of the paper. Although the Itô interpretation is used here, we could also consider a Stratonovich interpretation by modifying the particular details of $\mathbb{L}$, as this description still produces a differential operator in the forward equation.

2.3 Jump component

Let the jump process be an inhomogeneous (state-dependent) Poisson process with intensity (rate) $\lambda(X_{t})$. Also necessary is a description of what occurs at the jumps, sometimes referred to as the reset map [15]. The behavior of the reset map is characterized by the jump operator $\mathbb{J}$, which is a probability density flux, ensuring that (2) indeed describes the evolution of a probability density. We now describe a few possible choices of this operator, letting $\tau$ denote the jump time.

2.3.1 Jump operator examples

Constant jump size. At the jump times, the process increments by a fixed amount $\Delta$, so

[TABLE]

Then, the corresponding jump operator $\mathbb{J}$ is a shift by that fixed quantity

[TABLE] 2. 2.

Reset. Fundamentally different than the previous example, rather than jumping by a fixed displacement, the process resets to a specific position, $\eta$, so that

[TABLE]

The corresponding $\mathbb{J}$ operator is

[TABLE]

where $\delta(x)$ is the Dirac delta function. The integral scaling is necessary to preserve probability fluxes: $\int_{-\infty}^{\infty}\mathbb{J}p\,\text{d}x=\int_{-\infty}^{\infty}p\,\text{d}x$. 3. 3.

Random jump size. A generalization of the first example, the particle can jump a random displacement $\Delta$, where the size of the jump is described by the probability density $\Delta\sim\mu(\Delta,x)$, which can also be state-dependent. The $\mathbb{J}$ operator is

[TABLE] 4. 4.

Other examples. Although not directly considered in this paper, maps of the form

[TABLE]

have been studied elsewhere [15] and can be formulated in our framework in the described manner. A generalization of this case could also be made to allow for jumps to a random location, but this does not seem to provide additional interesting structure.

While diffusion with resetting fits within the proposed framework, this topic has a long, rich history of study [23, 45]. Included within this body of literature are investigations of: time-dependent resetting ($t$-dependent $\lambda$) [48], optimal resetting [22], resetting in bounded domains [9], resets with drift [46], characterizations of the stationary behavior [43], and path-integral formulations [50]. We do not seek to propose our framework as an alternative to this rich body of literature, but merely identify that stochastic resetting fits within a family of models that otherwise lack this same level of attention and can be studied with the tools presented here.

3 Results

In this section, we present the general theoretical results of the paper regarding the distributions of jump locations and interjump times.

3.1 Jump distributions

We begin by providing intuition for the construction of the sequence of jump locations. Rather than studying the full process including diffusion and jumps, a more convenient process to study is a survival formulation between jumps, which we denote the absorbing process,

[TABLE]

Note that we distinguish between this and the full process $\tilde{p}(x,t)$, which we refer to henceforth as the reinjected process. While the absorbing process does not capture all the behavior of the full re-injected process (as it does not contain any information about $\mathbb{J}$), the absorbing process fixes the distribution (in both space and time) of jump events and consequently serves more fruitful in disentangling the jump component from the full process.

Theorem 1**.**

The densities of the jump location $p_{\ell}(x)$ and jump time $p_{\tau}(t)$ of (11) are described by

[TABLE]

and

[TABLE]

Proof.

Define the survival probability $s(t)$ to be the probability that the jump has not occurred by time $\tau$

[TABLE]

which means the first jump time density $p_{\tau}(t)$ is then

[TABLE]

Then, using (11) and the fact that $p\to 0$ as $x\to\pm\infty$, we obtain

[TABLE]

Due to the lack of spatial flux in $q(x,t)$, the density of exit locations is $q(x,\infty)$, which can be obtained by integrating over all possible jump times

[TABLE]

by noting that $q(x,0)\equiv 0$. For this proof (and the remainder of the results), we assume that $\lambda,\mathbb{L},\mathbb{J}$ are chosen with sufficient regularity such that there are no explosions (i.e. the intensity stays finite. For more details on these conditions, see [3]. ∎

Although this provides useful information about the distributions (in time and space) at which a particular jump occurs, we are interested in studying the distributions of all jumps. This requires including the information embedded in $\mathbb{J}$ about how to re-inject the particle. triggered. A similar idea of studying the system between jumps was proposed in [47], in which the authors study the record statistics of the so-called Sisyphus random walk. However, the Sisyphus random walk is a discrete random walk on a lattice, rather than a continuous diffusion as discussed here, meaning the results are not directly applicable but strongly paralleled.

3.2 Jump location sequential mapping

We ystudy the interjump dynamics by noting that between jumps, the full process (2) can be described by the absorbing process (11). Define $t_{i}$ to be the $i$th jump time and let $X_{i}\coloneqq X_{t_{i}}$ be the $i$th jump location. Let $p_{i}(x)$ be the probability density of $X_{t_{i}}$ formally defined to be,

[TABLE]

Theorem 2**.**

The distribution of the $i$th jump location, $p_{i}(x)$ is described by the following iterative relation

[TABLE]

where $p_{0}$ be some known starting distribution of the process.

This construction follows from directly from Theorem 1. Effectively, $\hat{p}_{i}$ serves as intermediate quantity tracking the distribution of all possible jump locations (and times) for that particular iterate. That is, $\hat{p}_{i}$ is the survival density (in both time and space), tracking the process between jumps, so

[TABLE]

Then, to start the next iterate, the distribution of jump locations must be modified by the jump procedure characterized by $\mathbb{J}$. Although tracking the jump locations rather than the jump times seems counterintuitive at first, no natural analogous formulation is apparent.

Define the more convenient quantity to study

[TABLE]

If $\lambda(x)=0$ for some $x$, then necessarily $p_{i}(x)=0$ since $\lambda(x)=0$ implies no jump can occur at this location, hence this quotient causes no difficulties.

Theorem 3**.**

The description (19) is equivalent to

[TABLE]

or more explicitly

[TABLE]

Proof.

Integrating both sides of (19) with respect to $t$ and noting that $\hat{p}_{i}(x,\infty)=0$ we are left with

[TABLE]

The linear operator $\mathbb{L}$ is a differential operator in $x$ and consequently commutes with the time integral. Using the initial condition and $p_{i}/\lambda=\int_{0}^{\infty}\hat{p}_{i}\,\text{d}t$, we obtain the desired result. ∎

The map (23) provides an explicit description for the sequence

[TABLE]

where, we assume $u_{\star}$ is the fixed point of of the map (23) and exists. Although we seek to not dwell on this aspect, a natural question to ask is under what conditions does $u_{\star}$ exist? Equivalently, when does the the process reach stationarity? In B, we provide a brief commentary about how the relationship (23) can be thought of as an iterated linear non-negative integral operator which can be studied as such. Results from theoretical ecology literature are then cited which provide a heuristic analysis of conditions for convergence.

We proceed assuming that (23) has a fixed point, which must be of the following form.

Theorem 4**.**

The stationary distribution of the jump locations $p_{\star}(x)$ is described by

[TABLE]

where $u_{\star}\coloneqq p_{\star}/\lambda$.

Corollary 1**.**

The stationary distribution of the jump locations $p_{\star}$ is the same as the stationary distribution of the full process if and only if $\lambda$ is constant.

Proof.

This is an immediate consequence of taking the full process to be in stationarity, so $\text{d}\hat{p}/\text{d}t=0$ in (2), which satisfies

[TABLE]

This is exactly the same relationship as (26). Recalling that $\lambda u_{\star}=p_{\star}$ and that both $p_{\star}$ and $\hat{p}_{s}$ are probability densities, the only way that $p_{\star}=\hat{p}_{s}$ is if $u_{\star}$ is a rescaling of $p_{\star}$ or, equivalently, $\lambda$ is constant. ∎

Elaborating a bit on this relationship: while (27) and (26) may appear to produce the same result, the solutions to each require different scalings. Since $\hat{p}_{s}$ is a probability density, it must be that

[TABLE]

However, $u_{\star}\lambda=p^{\star}$, where $p_{\star}$ is a probability density, so this means that

[TABLE]

It is worth pointing out that [14] derives similar results relating the jump locations and jump times, however, the justification used there depends on the trajectories being deterministic between jumps, i.e. without diffusion.

3.3 Moments

Define $\tau_{i}$ to be the $i$th interjump time, $\tau_{i}\coloneqq t_{i}-t_{i-1}$. Also, let $\tau_{\star}$ be the interjump time for the stationary process.

Theorem 5**.**

The mean interjump time $\tau_{i}$ can be recovered from the distribution of the $i$th jump location and is

[TABLE]

Also, the mean stationary interjump time can be computed from the stationary distribution of the jump locations $p_{\star}$

[TABLE]

Proof.

We integrate both sides of (19) with respect to $x$, again noting that $\hat{p}_{i}\to 0$ as $x\to\pm\infty$, resulting in

[TABLE]

However, from Theorem 1, we see that the right-hand side is exactly the distribution of the interjump time $p_{\tau_{i}}$, so we have

[TABLE]

Taking the mean on both sides with respect to $\tau_{i}$,

[TABLE]

after integrating by parts and noting that $\int_{0}^{\infty}\hat{p}_{i}\,\text{d}t=p_{i}/\lambda$, we get the desired result. ∎

While the first moment (mean) of the interjump time can be computed directly with a quadrature of the presumed known jump location $p_{i}$, the higher order moments are less straightforward. A describes how, in theory, knowledge of the $p_{i}$ can be used to extract higher order moments of $\tau_{i}$. Solving for higher order moments using this approach requires solving a hierarchy of differential equations, which in practice, may not be so feasible. However, a more practical numerical approach may be to solve for $p_{i}$ and then run the absorbing (11) process to extract explicitly $p_{\tau_{i}}$ using the relationships from Theorem 1.

4 Examples

4.1 Neuronal integrate-and-fire

The behavior of individual neurons can roughly be thought of as an “integration” process, which builds up voltage, and then “fires,” which describes the release of this voltage into an action potential [28], the full process of which is known as integrate-and-fire. The most classical version of this model involves a deterministic buildup of voltage until a fixed firing voltage threshold is reached.

Incorporating noise in various ways into this class of models has a rich body of literature [51, 49]. Some approaches include a so-called “leaky” neuron, one that has diffusive noise in the integration phase [19, 16, 25], whereas others regard the threshold itself as stochastic [6, 13, 42, 19]. Our work provides a natural framework to consider both sources of noise and their effect.

We propose a simple model of neuronal integrate-and-fire. While this model is greatly simplified from the actual physiology, we take this approach to show that a minimal model, stripped of considerable details, is still able to produce interesting emergent behavior. In previous models with a stochastic threshold, there is some inherent threshold, say, $v_{0}$, where the firing rate $\lambda(v)$ is taken to be a Gaussian centered around $v_{0}$ [13, 6]. The justification for this that models such as Hodgkin-Huxley [35], which account for more fine-grained detail, predict a distinct firing threshold $v_{0}$. In our model, we find that this sharp threshold $v_{0}$ can also arise from stochasticity alone, even when no inherent threshold is defined in the model.

Let $V_{t}$ denote the voltage at time $t$. Then, in the context of this framework, we take the forms

[TABLE]

In words: we take the buildup to be at a constant rate $\alpha$ and leaky with constant diffusivity. While it is typically more common to take evolution of the form $\dot{v}=-a+bv$, this defines an inherent voltage $\tilde{v}=b/a$, which we deliberately omit. The result of making this change is also qualitatively negligible. For the jump component, we take $\mathbb{J}$ to be the reset to zero case, a full fire. The rate at which the firing occurs is assumed to be monotonically increasing, where again, we emphasize the fundamental difference between this form and previously considered, where $\lambda$ is a Gaussian centered around some voltage $v_{0}$. In our model, voltage simply builds up, and as it builds up, the system is more likely to fire. There is no prescribed firing threshold as in classical work and its extensions. The justification for this form can be thought of as the underlying mechanism for firing: a discrete birth-death process (also with voltage dependent rates) of ion channels. At some fixed voltage $v$, the birth-death process is more likely to flip entirely “on” and cause the neuron to fire. The relationship between models with explicit ion channels and integrate-and-fire has been explored previously [17].

An example simulation of the model can be seen in Figure 2(a) with parameter values $\alpha=2,D=.1,\gamma=.5,\beta=2$. In the simulation, voltage builds up in a noisy (leaky) manner, and then the jump process fires, resetting back to $v=0$, typical of an integrate-and-fire model. However, in this case, we emphasize that the firing voltage is inherently random and therefore of interest to study. Then, in stationarity, the fundamental relation in this work (26) becomes (noting $x\to v$)

[TABLE]

From (29), the integral scaling on the $\delta$ function is equal to $1$,

[TABLE]

This can be computed in a similar manner to a Green’s function, noting that for $v\neq 0$, we have

[TABLE]

and integrating (37) from $(-\varepsilon,\varepsilon)$ and using the fact that $u_{\star}$ must be continuous, we get the matching condition

[TABLE]

Solving (38), we get that our solution is of the form

[TABLE]

where $I,K$ are modified Bessel functions [2] such that $I\to 0$ as $v\to-\infty$ and $K\to 0$ as $v\to\infty$ (and are linearly independent). The jump condition (39) then becomes

[TABLE]

However, we also need to impose continuity of $u_{\star}$, so we also have the requirement

[TABLE]

The conditions (41),(42) provide us two equations for two unknowns, which yield

[TABLE]

Thus far, our solution is defined only up to a constant, but we know that $u_{\star}$ must satisfy the scaling (29), which our choice of $c_{1},c_{2}$ serendipitously already satisfy. By computing $u_{\star}$ explicitly, we can use (21) to immediately obtain the distribution of jump locations $p_{\star}$. Finally, from Corollary 1, we also have the stationary density of the full process, $\hat{p}_{s}$, which is a rescaling of $q_{\star}$ such that (28) is satisfied.

The stationary density $\hat{p}_{s}$ of the full process and the stationary jump distribution $p_{\star}$ are shown in Figure 2(b,c). From this, we can see that $\hat{p}_{s}$, the stationary density of the full process provides no interesting information. It is relatively uniform through some range of voltages. If a PDE approach were taken to provide information about the full process, this is all that would be available. However, the jump locations $p_{\star}$, are noteworthy. Despite prescribing no explicit voltage threshold for firing, the state-dependent nature of the firing rate yields an effective threshold ($v_{0}\approx 1$ for these parameter values) and persists through a wide range range of parameters. Hence, this is another, novel justification for models that take a deterministic firing threshold, as our minimal integrate-and-fire model produces this feature as a product of stochasticity, which is made apparent in the lens of the proposed framework.

4.2 Transport by a molecular motor

Consider a single molecular motor (e.g. kinesin [31]) attached to a cargo as in Figure 3(a). Molecular motors produce transport by taking discrete steps along a track, exerting a force on a cargo. A notable feature of these steps is that the rate at which they are taken is well-established to be force-dependent [8]. That is, motors typically “slow down” as force is exerted on them. As this discrete stepping process occurs, another source of noise is the diffusion of the cargo, which is also instantaneously changing the force applied to the motor, and consequently, the rate at which it steps. In our proposed simple motor of transport, we provide a preliminary investigation into this question.

We assume that the linker between the molecular motor and the cargo can be treated as a Hookean spring, where $x$ denotes the distance stretched (from rest), as depicted in Figure 3(a). Then, the diffusion component of this process is an Ornstein-Uhlenbeck process. We take the motor to take fixed step sizes, $\Delta$ and the rate at which motors step to be Gaussian $\lambda(x)=\sqrt{\frac{\alpha^{2}\beta}{2\pi}}\exp\{-\beta x^{2}/2\}$. Although there is evidence that the direction of the force is significant for motor stepping, this often deals with the proportion of forward to back steps of the motor. Because we assume forward steps here, we assume, for simplicity, any force decreases the motor stepping rate.

The resulting stationary relationship (26) becomes

[TABLE]

There does not appear to be a fruitful approach to studying (44) analytically. However, we solve can it numerically using a Crank-Nicolson upwind scheme. From this, $u_{\star}$ and consequently $p_{\star}$ is obtained. The first moment of the interjump time $\langle\tau_{\star}\rangle$ is immediate after a quadrature (31).

We plot the mean value of $x$ when a jump occurs, denoted $\langle j\rangle$ and the mean time between jumps, $\langle\tau\rangle$ as a function of the cargo diffusion coefficient $D$, for varying values of the linker strength, $a$. These results can be seen in Figure 3(b,c) for parameter values $D=.1,\alpha=1,\beta=5$. The strength of diffusion appears to produce a curious effect on the other noise source. While the mean jump location, $\langle j\rangle$, remains effectively unchanged for varying values of linker strength, the mean interjump time exhibits very different behavior for different values of $a$. Specifically, we see that if the linker is weak, the interjump time can have a non-monotonic dependence on the diffusion coefficient. That is, the interaction between these two sources of noise (which can be studied as a result of this work) is non-trivial.

While this model is far simpler than the complex physiological it describes, the overall result (non-monotonic dependence on cargo diffusion) may be a feature that persists in a more detailed motor model (e.g. multi-motor), as other non-monotonic dependences on cargo noise have been seen [44]. Although the quantity $\langle j\rangle$ (the average distance a motor is stretched when stepping) is perhaps not the most useful in extracting information about more complex transport, it seems feasible that it it could be utilized in a reward-renewal framework akin to others utilized to study motor transport such as [40, 32, 33]. Motors also often have an asymmetric (state-dependent) unbinding rate, thought to be significant [53] and could perhaps provide a future application for this framework.

4.3 Independence of $\lambda(x)$

Lastly, we provide a short example (with no physical motivation in mind) that shows that our framework can be used with other tools (here, Fourier analysis) to glean interesting features of certain families of models. For this example, we consider a constant drift and diffusion with random jump size distribution $\mu(\Delta)$. Then, in stationarity, (26) becomes

[TABLE]

One quantity of interest is the first moment of jump locations, so

[TABLE]

Although $\lambda(x)$ is arbitrary, surprisingly the Fourier transform can be used to gain some insight. Define the Fourier transform of $u_{\star}$ to be

[TABLE]

Also define the transformed quantities

[TABLE]

from which, we note that

[TABLE]

Taking the transform of (45),

[TABLE]

which, evaluated at $k=0$ yields

[TABLE]

However, we know $\mu(\Delta)$ and $u_{\star}\lambda=p_{\star}$ are both probability densities and consequently $M(0)=1$ and $\Lambda(0)=1$. Thus, (51) is trivially satisfied. This is not useful on its own, but we can couple the higher order moments by taking a $k$ derivative of (50) to yield

[TABLE]

which, at $k=0$ and using $M(0)=1,\Lambda(0)=1$ yields

[TABLE]

which says that necessarily

[TABLE]

However, recall from (30) that $U(0)=\langle\tau_{\star}\rangle$, the mean interjump time, so

[TABLE]

This result is interesting for two reasons. For one, we have the mean interjump time is completely independent of the choice of $\lambda(x)$, aside from the requirement that the system reaches stationarity. Secondly, the sign of (55) narrows down the requirements for the process to reach stationarity. Since $\langle\tau_{\star}\rangle$ is an interjump time, it must be non-negative, meaning this quantity only exists if $a,\langle\mu\rangle$ differ in sign. This is intuitive as the drift and jump process must oppose each other to have a chance of reaching stationarity. We acknowledge that the linearity of this example (which allows us to perform the explicit Fourier analysis) also makes the result possible to derive through other methods (e.g. taking expectations of the SDE). However, we hope this analysis conveys that our framework provides a different lens (in junction with elementary tools) to study and possibly reveal behavior of a system.

5 Discussion & Conclusion

In this work, we present a general framework for studying jump-diffusion systems with state-dependent jump rates, a class of models that were previously difficult or impossible to study otherwise due to the interaction between the two sources of noise. The formulation is flexible enough to accommodate a variety of behaviors, providing relevance to a wide range of applications. The particular objects of study in this work are the distributions of the jump locations: the values of $X_{t}$ evaluated at the jump times $t_{i}$.

We reformulate the full process into a survival description between jumps. Using this reformulation, the explicit distribution of jump locations can be extracted in a manner that depends on the previous jump location. This relation can be articulated as an iterated map, producing a sequence of jump locations. With explicit knowledge of the distribution of jump locations, statistics of the interjump times can be computed in more generality than previous work [14]. Taking this map to its limit, or equivalently, assuming the process reaches stationarity, the stationary distribution of the full process and the jump locations are shown to be closely related. An immediate consequence of this relation is that these two quantities differ if and only if the jump intensity is state dependent. That is, for the class of models studied in this work, the jump locations are distinct from the stationary density. Consequently, if the jump locations are of interest, the results of this work illustrate their relation to observations of the full process.

Although not the focus of this work, we provide a brief discussion of possible conditions on convergence to stationarity in the novel lens of this framework. This discussion utilizes tools from theoretical ecology motivated by the observation that the iterated map formulation is effectively a non-negative integral linear operator. This is in contrast to previous, more probabilistic approaches to studying finite time blowup of these models or similar [3, 55, 15, 10].

Finally, we provide three example applications of the framework. The first application, to a simple model of stochastic neuronal integrate-and-fire, is used to show a scenario where the jump locations themselves may be of interest. Our minimal model contains no inherent threshold voltage, but, as a product of stochasticity of the firing threshold, forms a robust, sharp peak of firing locations. This provides additional and novel support for the use of a deterministic threshold as in classical integrate-and-fire and its extensions. The second model, one of intracellular transport by a molecular motor is used to demonstrate how the two sources of noise (diffusion and jumps) interact. We ultimately find that the mean stepping rate of the motor can have a non-monotonic dependence on the strength of cargo diffusion, a consequence of the interaction between the two noise sources. Finally, a third example calculation shows that certain classes of models (in this case, constant drift, random jump size) can have curious behaviors. Specifically, we find that the mean interjump time for be independent of the choice of $\lambda(x)$, so long as the process does indeed reach stationarity.

This work proposes a preliminary framework for studying a wide class of problems. Although not explicitly studied here, we suspect that with more exotic diffusions (e.g. colored noise, which fits within the scope of these results), the interactions between the feedback between the noise sources can produce even more interesting behavior and will be explored in the future. Additionally, a more precise exploration is necessary of how an iterated map framework such as this can be used to provide conditions of stationarity.

6 Acknowledgments

This research was partially supported by NSF grant DMS 1515130 and DMS-RTG 1148230.

Appendix A Higher order interjump time moments

In theory, higher order moments of the interjump times can be computed with knowledge of the sequence of jump distributions $p_{i}$ as was discussed with the first moment in (30).

To illustrate this, consider (33) and take the second moment with respect to $\tau_{i}$ on both sides, resulting in

[TABLE]

Again, noting that $u_{i}=p_{i}/\lambda=\int_{0}^{\infty}\hat{p}_{i}\,\text{d}t$, after integrating by parts, we get

[TABLE]

which looks similar to (30), however, $v_{i}$ is defined by the relation

[TABLE]

Thus, computation of the second moment requires solving the differential relationship (58), rather than just a quadrature as in the first moment. Continuing to higher order moments using same approach yields a further hierarchy coupled in a differential manner. Interestingly, the differential operator (left hand side) of (58) is exactly $\mathbb{T}$, the same as (23), meaning the Green’s function (effectively $\mathbb{T}^{-1}$, discussed more in B) could be used to compute these quantities.

Appendix B Spectral properties of iterated map

Note that (23) involves the inverse of a differential linear operator $\mathbb{T}$, which is exactly determined by its corresponding Green’s function. Let $G(x,\xi)$ be the corresponding Green’s function to $\mathbb{T}$, meaning that

[TABLE]

Then, (23) can be rewritten as

[TABLE]

After a change of variables, define $\tilde{G}$ by

[TABLE]

For example, for fixed jump sizes: $\mathbb{J}p=p(x-\Delta)$, then

[TABLE]

Abbreviate this linear non-negative integral operator

[TABLE]

Now, iterations of our map correspond to iterating the integral operator $\mathbb{A}$ with kernel $\tilde{G}$. This formulation is exactly that of the so-called integral projection models (IPM) in theoretical ecology[21]. Stability results from the IPM literature can be used to understand the convergence of our iterative procedure.

The main result comes from [20] and is effectively a statement of the Krein-Rutman theorem [39, 38], the infinite dimensional analog of the Perron-Frobenius theorem. [34, 7]

Although $L^{1}$ seems like the natural function space to study the spectral properties of these operators, since they must preserve probability, $L^{2}$ turns out to be far more accessible due to issues establishing compactness in $L^{1}$. In Appendix C of [21], the authors provide a more thorough commentary on these complications. we cite the main theorem which establishes the existence of a dominant eigenvalue for $\mathbb{A}$.

Theorem B.1** (Easterling 1998, [20]).**

Suppose that $\tilde{G}\in L^{2}$ and is non-negative. If there exists an $\alpha>0,\beta>0,u_{0}$ such that

[TABLE]

for all $x,\xi$, then the integral operator $\mathbb{A}$ has a dominant eigenvalue with associated eigenfunction.

This establishes conditions for the existence of a dominant eigenvalue and eigenvector, which establish the long-term behavior.

Theorem B.2** (Easterling 1998, [20]).**

Assuming that $\mathbb{A}$ satisfies the previous condition, then the stationary distribution $u_{\star}$ is given by the eigenfunction $\phi_{1}$ associated with the dominant eigenvalue $\lambda_{1}$,

[TABLE]

where $\kappa$ is a scaling parameter.

The previous theorem suggests that the sequence $u_{i}$ converges if and only if $|\lambda_{1}|=1$. Thus, this summary provides a rough heuristic, but possibly novel angle for determining whether the sequence $u_{i}$ (and the full process) approach stationarity, although we reemphasize that this is not our focus.