2
$\begingroup$

I need to find the distribution of the random variable $T_{peak}$ where $T_{peak}$ represents the peak time of the first hitting time process.

Detailed Explanation of the System: There are $N^{Tx}$ amount of emitted molecules from a specific point in 3D environment. Molecules diffuse in the environment according to the followings:

$$ r[t] = r[t-1] + (\Delta r_1, \Delta r_2, \Delta r_3)$$ $$ \Delta r_i \sim \mathcal{N}(0,\, 2D\Delta t)$$ where $r[t]$, $r_i$, $D$, and $\Delta t$ are the location vector at time $t$, $i$-th component of the location vector, diffusion coefficient, and the time step, respectively.

If there is an absorbing spherical trap at a distance $d$, the mean number of arriving/hitting molecules until time $t$ is:

$$ E[N^{absorb}(t)] = N^{Tx} \frac{r_{trap}}{d+r_{trap}} \, \text{erfc} \left( \frac{d}{\sqrt{4Dt}} \right) = N^{Tx} \frac{r_{trap}}{d+r_{trap}} \, 2\Phi \left( \frac{-d}{\sqrt{2Dt}} \right)$$ where $r_{trap}$ is the radius of the absorbing spherical trap.

When you plot $N^{absorb}(t)$ in small intervals, you get something like the following figure (Scaled Inverse Gaussian distribution) enter image description here

And the expected value of the peak time of hitting time histogram is $$ E[T_{peak}] = \frac{d^2}{6D}$$

When simulating this diffusion process and the focusing on the absorbing times, $T_{peak}$ differs from simulation instance to instance, $T_{peak}$ is a random variable and I need to find the distribution of $T_{peak}$.

::Edit/New Text Start

If we consider time bins $\{t_0, t_1, ..., t_k, ...\}$ and denote the number of absorbed molecules in that bin as $N(t_k)=N^{absorb}(t_k^-)-N^{absorb}(t_k^+)$ where $t_k^-$ and $t_k^+$ are the left and the right end of the $k$-th time bin, we can define $T_{peak}$ as follows: $$T_{peak} = \arg\max\limits_{t_k} N(t_k)$$ with this definition, what is the distribution of $T_{peak}$. Is it Poisson or Binomial or Gaussian or Inverse Gaussian? Does it have well known structure and name?

Example Topology Figure: enter image description here

::Edit/New Text End

P.S. These absorbed molecules are considered as the received signal in molecular communications and the peak time distribution of the received signal is important for many applications.

$\endgroup$
  • $\begingroup$ - r:location vector has three components (r_1, r_2, r_3) you may consider that as (x,y,z) $\endgroup$ – fermat4214 Jan 30 at 13:47
  • $\begingroup$ - r_trap is the radius of the fully absorbing sphere so the distance between the emission point and the center of the absorbing sphere is d+r_trap $\endgroup$ – fermat4214 Jan 30 at 13:48
  • $\begingroup$ That is for single molecule emission, I will and $N^{Tx}$ number of emitted molecules $\endgroup$ – fermat4214 Jan 30 at 14:44
  • $\begingroup$ I added $N^{Tx}$ for convenience, thx $\endgroup$ – fermat4214 Jan 30 at 14:47
  • $\begingroup$ Let us continue this discussion in chat. $\endgroup$ – fermat4214 Jan 30 at 15:00
2
+50
$\begingroup$

Description as multinomial distribution

You could view the distribution of the particles in the different time bins as a multinomial distribution where the mean amount of particles $\lambda_k$ in a particular time bin $t_k$ is $N(t_{k+1})-N(t_{k})$, or approximating this with the derivative of $N(t)$

$$\lambda_k = N^\prime(t) \Delta t = \frac{r_{trap}}{d+r_{trap}} \frac{n}{\sqrt{2Dt_k^3}} \phi\left(\frac{-d}{\sqrt{2Dt_k}}\right) \Delta t$$

where $n$ is the total number of particles.


Approximations

When the number of bins is large, when there is less correlation between pairs of bins, then you can approximate the multinomial distribution by individual binomial distributions.

When the number of particles is large then these binomial distributions can be approximated by poisson distributions or normal distributions.


Probability distribution for time bin with maximum number of molecules

With the distribution of number of molecules $n_k$ in the $k$-th time bin described by:

$$n_k \sim Pois(\lambda_k)$$

You can compute the probability for a particular bin to contain the extreme value as

$$P(\text{max at $t_k$}) = \sum_{i=0}^\infty \underbrace{f(i,\lambda_k)}_{\substack{\text{probability bin $k$} \\ \text{ contains $i$ molecules}}} \prod_{j \neq k} \underbrace{ F(i,\lambda_j) }_{\substack{\text{probability bin $j$} \\ \text{ contains at most $i$ molecules$$}}} $$

Where $f$ and $F$ are the density and cumulative distribution function of the Poisson distribution.


This is a bit of an horrible expression, but it can already be computed. I will see at a later time If I can make it look easier.

Computation example:

The image below shows the outcome of a simulation along with the theoretic distribution for the peak time (which seems to work well)

computation

It resembles a curve that might be approximated with a normal curve (although it is not entirely the same).

qqplot

# model parameters
dt <- 0.001
t <- seq(dt,0.3,dt)
n=700*1000
D = 1
d = 0.5

# model
ft <- n*d/sqrt(2*D*t^3)*dnorm(d/sqrt(2*D*t),0,1)
fmids <- n*d/sqrt(2*D*(t+dt/2)^3)*dnorm(d/sqrt(2*D*(t+dt/2)),0,1)
plot(t,ft*dt,type="l",lwd=1.5,lty=2)

# simulation
#    
# simulation by drawing from uniform distribution
# and converting to time by using quantile function of normal distribution
ps <- runif(n,0,1)                
ts <- 2*pnorm(-d/sqrt(2*D*t))     
sumn <- sapply(ts, FUN = function(tb) sum(ps < tb))
lines(t[-length(sumn)],sumn[-1]-sumn[-length(sumn)],col=4)

# compute t peak distribution
pmax = rep(0,length(t))         # probability of max
for (i in 1:3000) {
  Fj <- ppois(i,fmids*dt,log.p = TRUE)
  ss <- sum(Fj)
  for (j in 1:length(t)) {
    # doing the product by taking sum of log values
    pmax[j] <- pmax[j]+exp(dpois(i,fmids[j]*dt,log = TRUE) + ss - Fj[j])
  }
}

# compute whether it sums to 1
# you will see a bit higher because a maximum value may occur at multiple times
sum(pmax)

#multiple simulations
nsim <- 10000
histt <- c(0) #container with simulated maximums
pb <- txtProgressBar(title = "progress bar", min = 0,
                     max = nsim, style=3)
for (i in 1:nsim) {
  ps <- runif(n,0,1)
  ts <- 2*pnorm(-d/sqrt(2*D*t))
  sumn <- sapply(ts, FUN = function(tb) sum(ps < tb))
  histt <- c(histt,which(sumn[-1]-sumn[-length(sumn)] == 
                            max(sumn[-1]-sumn[-length(sumn)])))
  setTxtProgressBar(pb, i)
}
close(pb)
histt <- histt[-1]  # strip leading 0

h <- hist(histt*dt,breaks = t,xlim=c(0.025,0.06),main="peak time",
          xlab = "time")
lines(h$mids,pmax[-1]*nsim,col=2)
points(h$mids,pmax[-1]*nsim,col=2)

qqnorm(histt, pch=21,cex=0.5,col=1,bg=1)
qqline(histt)
$\endgroup$
  • $\begingroup$ Let me again clarify two points. First, r:location vector has three components (r_1, r_2, r_3) you may consider that as (x,y,z). Second, r_trap is the radius of the fully absorbing sphere so the distance between the emission point and the center of the absorbing sphere is d+r_trap. Now according to your notation let me rephrase my question: $T_{peak} = \arg\max\limits_{t_k} N(t_k)$ with this definition, what is the distribution of $T_{peak}$. Is it Poisson or Binomial or Gaussian or Inverse Gaussian? Does it have well known structure and name? $\endgroup$ – fermat4214 Jan 30 at 14:00
  • $\begingroup$ Molecules are emitted from a point source, you may assume any point in space. There is a spherical absorbing trap at distance $d$ and the radius of the spherical absorber is $r_{trap}$. Please note that the distance between emission point and the center of the absorbing sphere is $d+r_{trap}$ $\endgroup$ – fermat4214 Jan 30 at 14:14
  • $\begingroup$ I have added an example topology figure to the question body. $\endgroup$ – fermat4214 Jan 30 at 14:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.