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)
Code: Select all
bkgdegredation=rxd.Reaction(tfact>tfactprod, 1e-3)
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