#!/usr/bin/env python
# coding: utf-8

# # Reversible (Diffusion-limited)
# 
# This is for an integrated test of E-Cell4. Here, we test a simple reversible association/dissociation model in volume.

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
from ecell4.prelude import *


# Parameters are given as follows. `D`, `radius`, `N_A`, `U`, and `ka_factor` mean a diffusion constant, a radius of molecules, an initial number of molecules of `A` and `B`, a ratio of dissociated form of `A` at the steady state, and a ratio between an intrinsic association rate and collision rate defined as `ka` and`kD` below, respectively. Dimensions of length and time are assumed to be micro-meter and second.

# In[2]:


D = 1
radius = 0.005
N_A = 60
U = 0.5
ka_factor = 10  # 10 is for diffusion-limited


# In[3]:


N = 20  # a number of samples


# Calculating optimal reaction rates. `ka` and `kd` are intrinsic, `kon` and `koff` are effective reaction rates.

# In[4]:


import numpy
kD = 4 * numpy.pi * (radius * 2) * (D * 2)
ka = kD * ka_factor
kd = ka * N_A * U * U / (1 - U)
kon = ka * kD / (ka + kD)
koff = kd * kon / ka


# Start with no `C` molecules, and simulate 3 seconds.

# In[5]:


y0 = {'A': N_A, 'B': N_A}
duration = 0.35
opt_kwargs = {'legend': True}


# Make a model with effective rates. This model is for macroscopic simulation algorithms.

# In[6]:


with species_attributes():
    A | B | C | {'radius': radius, 'D': D}

with reaction_rules():
    A + B == C | (kon, koff)

m = get_model()


# Save a result with `ode` as `obs`, and plot it:

# In[7]:


ret1 = run_simulation(duration, y0=y0, model=m)
ret1.plot(**opt_kwargs)


# Make a model with intrinsic rates. This model is for microscopic (particle) simulation algorithms.

# In[8]:


with species_attributes():
    A | B | C | {'radius': radius, 'D': D}

with reaction_rules():
    A + B == C | (ka, kd)

m = get_model()


# Simulating with `spatiocyte`. `voxel_radius` is given as `radius`. Use `alpha` enough less than `1.0` for a diffusion-limited case (Bars represent standard error of the mean):

# In[9]:


# alpha = 0.03
ret2 = ensemble_simulations(duration, ndiv=20, y0=y0, model=m, solver=('spatiocyte', radius), repeat=N)
ret2.plot('o', ret1, '-', **opt_kwargs)


# Simulating with `egfrd`:

# In[10]:


ret2 = ensemble_simulations(duration, ndiv=20, y0=y0, model=m, solver=('egfrd', Integer3(4, 4, 4)), repeat=N)
ret2.plot('o', ret1, '-', **opt_kwargs)