RXD lower concentration limit?

Extending NEURON to handle reaction-diffusion problems.

Moderators: hines, wwlytton, ramcdougal

Post Reply
RBJ
Posts: 62
Joined: Sun Aug 02, 2015 4:28 am
Location: UK
Contact:

RXD lower concentration limit?

Post by RBJ »

Hello again, I m having good success with the RXD now, and hoping to incorporate in a range of models now. I’ve hit another hiccup however.
When I have a species that should be decreasing in concentration asymptotically towards zero it invariably equilibrates at somewhat above zero. For example, dropping from 2mM to 10uM.
So I read the post about setting absolute Rxd tolerance and added this line:

Code: Select all

tfact = rxd.Species(cytosol, d=0.5, initial=2,atolscale=1e-18)
And yet still with a simple degrading reaction... It never gets below about 10uM.

Code: Select all

bkgdegredation=rxd.Reaction(tfact>tfactprod, 1e-3)
So if I manually increase h.dt it gets me further down... But that seems I'll advised. If I just activate cvode I seem to lose stability completely.
Any tips please how I should approach this?
OK, will leave the above in place, but I realise it could, of course, be something weird I am doing so have put an example script of the phenomenon (incidently "<>" always throws a syntax error on my system, but that is unassociated with this issue)?

Code: Select all

import sys
print('Python version {}.{}'.format(sys.version_info.major,sys.version_info.minor))
import neuron
print ('Neuron Version',neuron.__version__)
from neuron import h, rxd
from matplotlib import pyplot
import numpy as np


h.load_file('stdrun.hoc')
cvode=h.CVode()
sec = h.Section(name='sec')
print("ATOL=",cvode.atol(1e-9))
#cvode.atol(1e-9)
sec.nseg = 11
sec.L=20
sec.diam=10
print ('cell length = {} um, max diameter {:.2f}um'.format(sec.L,sec.diam))
bigmax=0
rxd.set_solve_type(dimension=3)
cytosol = rxd.Region(h.allsec(),geometry=rxd.inside,dx=2)
tfact = rxd.Species(cytosol, d=0.5, initial=2.1,atolscale=1e-9)
tfactprod = rxd.Species(cytosol, d=0.5, initial=1e-9,atolscale=1e-9)
b = rxd.Species(cytosol, d=0.5, initial=1e-9,atolscale=1e-9)
print("ATOL=",cvode.atol())

bkgdegredation=rxd.Reaction(tfact,tfactprod + b, 1e-3)

h.finitialize(-60)  

#we will want an xy slice so where z=as close to zero as possible
minz=10000
for node in tfact.nodes:
    if 0<= node.z3d < minz:
        minz=node.z3d

#read these concentrations on demand
# I did it like this since data=tfact.nodes.value_to_grid() just gives me an array of NaNs
def gettimeslice():
    xt=[]
    yt=[]
    conct=[]
    for node in tfact.nodes:
        if (node.z3d==minz):
            conct.append(node.concentration)
            xt.append(node.x3d)
            yt.append(node.y3d)
    xt=np.array(xt).reshape(-1,1) 
    yt=np.array(yt).reshape(-1,1) 
    conct=np.array(conct,dtype=np.float64).reshape(-1,1) 
    
    return [xt,yt,conct]


from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

import pandas as pd

colors = [(0, 0, 1), (0, 1, 0), (1, 0, 0)]  # R -> G -> B
n_bins = [j for j in range(1, 500,10)]  # Discretizes the interpolation into bins
cmap_name = 'my_list'
cm = cm.colors.LinearSegmentedColormap.from_list(cmap_name, colors, N=len(n_bins))
switchONon=False
switchOFFon=False
switchCVODEon=False
ltime=0
inc=1

for i in range(1,20000,100):
    ltime=i*inc+ltime
    h.continuerun(ltime)
        
    #read the concentrations at z=0ish
    xm,ym,concm= gettimeslice() 

    # get the max on this loop, and a few other metrics for diagnosis
    maxconc=np.max(concm)
    minconc=np.min(concm)
    meanconc=np.mean(concm)
    #keep the biggest max ever in case missed it.
    if maxconc>bigmax:
        bigmax=maxconc
    print('T={:.2f}s, max={:.5f}mM, min={:.5f}mM, average={:.5f}mM, bigmax={:.5f} mM'.format(h.t,maxconc,minconc,meanconc,bigmax))
    
    #plot it nicely
    df1 = pd.DataFrame({'x': xm[:,0], 'y': ym[:,0], 'z': concm[:,0]})
    pyplot.close("all")
    fig = pyplot.figure()
    ax = Axes3D(fig)
    #ax.set_zlim(0, 0.25)
    
    surf = ax.plot_trisurf(df1.x, df1.y, df1.z, cmap=cm,vmin=0,
                            linewidth=0.1)
    cbar=fig.colorbar(surf, shrink=0.5, aspect=5)

Kind regards
Richard
ramcdougal
Posts: 267
Joined: Fri Nov 28, 2008 3:38 pm
Location: Yale School of Public Health

Re: RXD lower concentration limit?

Post by ramcdougal »

This is likely due to a previously under-appreciated detail in the numerical methods used for 3D simulation.

In particular, we use scipy.sparse.linalg.bicgstab with a tolerance of 1e-5 to solve a matrix problem.

This should have been connected to the atol used elsewhere*, but it is not. For now, you can choose a different tolerance via:

Code: Select all

MY_TOL = 1e-8
import scipy
from neuron import rxd
rxd.rxd._scipy_sparse_linalg_bicgstab = lambda a, b, M=None: scipy.sparse.linalg.bicgstab(a, b, tol=MY_TOL, M=M)
where MY_TOL is chosen as appropriate.

We expect that this issue will no longer be an issue in NEURON 7.7, as we have a version in development that gets better performance with a non-iterative solver.

Changes in Python syntax from Python 2 to Python 3 no longer allow using <> to separate two sides of a reaction. Instead, you can separate your left hand side and right hand side with commas:

Forward reaction:

Code: Select all

bkgdegredation=rxd.Reaction(tfact, tfactprod, 1e-3)
Bidirectional reaction:

Code: Select all

bkgdegredation=rxd.Reaction(tfact, tfactprod, 1e-3, 1e-5)
* Incidentally, cvode.atol and atolscale only affect variable step simulation, and you are using fixed step; either way, it would not help with the issue you are seeing.
RBJ
Posts: 62
Joined: Sun Aug 02, 2015 4:28 am
Location: UK
Contact:

Re: RXD lower concentration limit?

Post by RBJ »

Thank you yes that nails it!
Post Reply