PyEscape - A narrow escape simulation library

What is the narrow escape problem

Escape time solution for a unit sphere with N pores

\[\begin{aligned} f(\epsilon) &= \epsilon - \frac{\epsilon^2}{\pi} \log \epsilon + \frac{\epsilon^2}{\pi} \log 2 \\ \sigma &= \frac{N \epsilon^2}{4} \\ \kappa(\sigma) &= \frac{4\sigma}{\pi - 4 \sqrt{\sigma}} \\ \tau &= \frac{f(\epsilon)}{3D\kappa(\sigma)} + (15D)^{-1} \end{aligned} \]

\[\begin{aligned} D = \textrm{Diffusion coefficient} \end{aligned} \] \[\begin{aligned} \epsilon = \textrm{Pore radii} \end{aligned} \] \[\begin{aligned} N = \textrm{Number of pores} \end{aligned} \]

Assumes uniformly distributed pores

Narrow-escape problem for the unit sphere: Homogenization limit, optimal arrangements of large numbers of traps, and the N2 conjecture. Cheviakov et al. 2013

The simulation framework

The "escape" algorithm

  • We use time step sizes of 1e-8 seconds
  • Run 1000 simulations to get mean escape time
  • For one experiment of 1 second real time: 100000000000 (100 billion) movement calculations

The "escape" algorithm

        
def escape(r, cur_pos, delta, pore_locs, pore_size, check_func, max_steps=1e10, tol=1, dt=1e-8):
    pot_steps = fibonacci_spheres(max_steps, delta)
    idx = 0
    steps=0
    while idx < len(pot_steps):
        new_pos = cur_pos+pot_steps[idx]
        while (not check_func(new_pos, r)) and idx < max_steps: # while not in container (invalid movement)
            if passthrough_pore(new_pos, pore_locs, r=pore_size, tol=tol):
                    return (steps+1)*dt
            idx += 1
            new_pos = cur_pos+pot_steps[idx]
        cur_pos = new_pos
        steps+=1;idx+=1
    return -1

def check_func(p, r):
    return np.linalg.norm(p) < r

def passthrough_pore(p,p0, r=1, tol=1):
    return np.any(np.linalg.norm(p-p0, axis=1) < r*tol)
                    

1 https://github.com/AoifeHughes/NarrowEscapeSimulator/blob/master/PyEscape/escape_plan.py

More code examples

Multithreading and code reuse for multiple models:

        
                        from common_params import D, r, v, a, N, dt, MAX_steps, PS, ABCS

                        def esc(i):
                            p  = i[1]
                            pores = fibonacci_spheres(p, v)
                            np.random.seed()
                            res, loc = escape(D, v, a, pores, dt=dt, tol=1, max_steps=MAX_steps, with_exit_loc=True)
                            return res, loc
                        
                        if __name__ == "__main__":
                            idx = int(argv[1])
                            idy = 1
                            setups = {}
                            for p in PS:
                                for ABC in ABCS:
                                    setups[idy] = (p, ABC)
                                    idy += 1 
                            try:
                                p, ABC = setups[idx]
                            except:
                                print('Setup not found... exiting')
                                exit()
                            # this change here is to make comparable to the cubeoid experiments
                            p = len(points_on_cuboid(np.array(ABC), vol=v, npts=p))
                            args = [[i, p] for i in np.arange(0,N)]
                            with multiprocessing.Pool(processes=os.cpu_count()) as pool:
                                res = list(tqdm.tqdm(pool.imap(esc, args), total=N))
                        
                            res_t = [r[0] for r in res]
                            res_loc = [r[1] for r in res]
                            np.savetxt(f"./results/NEP_Sphere_{p}Pores_{N}Reps_{D}D_{v}V_{a}A_{dt}dt.csv", res_t)
                            np.savetxt(f"./exits/NEP_Sphere_{p}Pores_{N}Reps_{D}", res_loc)
                        

Turning into a command line app:

        
                            def main():
                                parser = argparse.ArgumentParser()
                                parser.add_argument('-D', help='Diffusion coefficient')
                                parser.add_argument('-v', help='Container volume')
                                parser.add_argument('-a', help='Escape pore area size')
                                parser.add_argument(
                                    '-s', help='Shape to escape (cube or sphere)', default='cube')
                                parser.add_argument('-p', help='Number of pores', default=1)
                                parser.add_argument(
                                    '-N', help='Number of simulations to run', default=1)
                                parser.add_argument(
                                    '-dt', help='Difference in time to use', default=3e-8)
                                parser.add_argument(
                                    '--cpu', help='Number of cores to use', default=1)
                                parser.add_argument(
                                    '-o', help='Output file', default="./results.csv")
                                args = parser.parse_args()
                                args = vars(args)
                                D, v, a, s, p, n, dt, cpu = float(args['D']), float(args['v']), float(
                                    args['a']), str(args['s']), int(args['p']), int(args['N']), float(args['dt']), int(args['cpu'])
                            
                                res, _ = run_simulations(D, v, a, s, p, n, dt, cpu)
                                print(np.mean(res))
                                np.savetxt(args['o'], np.array(res).flatten(), delimiter=',')
                                    
                        if __name__ == '__main__':
                            main()
                        

Turning into a command line app:

        
                        $ ./pyescape -D 100 -v 1, -a 0.1 -p 10 -N 100 --cpu 8 -o data.csv
                        

Analysis of data

To determine if the NEP solution was accurate in estimating escape time we build a series of random walk models

The NEP solutions show good agreement with our random walk models

With other containers and parameters the NEP solutions show robustness in estimating the mean escape time

Narrow escape solutions deviate with asymmetric shapes

Narrow escape solutions deviate with asymmetric shapes

Adjusting the NEP solution to use a larger volume, we find that a correction factor can be applied to cope with these asymmetric shapes

Narrow escape solutions deviate with asymmetric shapes

Adjusting the NEP solution to use a larger volume, we find that a correction factor can be applied to cope with these asymmetric shapes

Running these simulations we find that NEP solutions are in agreement with our models, but for higher number of pores this begins to diminish

Running these simulations we find that NEP solutions are in agreement with our models, but for higher number of pores this begins to diminish

Running these simulations we find that NEP solutions are in agreement with our models, but for higher number of pores this begins to diminish

Running these simulations we find that NEP solutions are in agreement with our models, but for higher number of pores this begins to diminish

Running these simulations we find that NEP solutions are in agreement with our models, but for higher number of pores this begins to diminish

Running these simulations we find that NEP solutions are in agreement with our models, but for higher number of pores this begins to diminish