I have programmed a model in both C++ and Python. This model has a noisy-input component, which I can replace with this C++:
double doubleRand() {
thread_local std::mt19937 generator(std::random_device{}());
std::normal_distribution<double> distribution(0.0, 1.0);
return distribution(generator);
}
Or this Python:
Inoise = (np.random.normal(0, 1) * knoise * np.sqrt(gNa * A))
IIon = ((iNa + iK + iL) * A) + Inoise #
# Compute change of voltage
v[i + 1] = (vT + ((-IIon + IStim) / C) * dt)[0]
The following is very strange: If I omit the noisy component (Inoise=0), then both models (C++ as well as Python) give exactly the same result. If I only introduce the noisy component (Istim=0), then both models give results (i.e. natural fluctuations that hardly differ from each other at 1000 runs). However, if I choose Istim=0.000001 and add noise, then the results differ by 30%. How is that possible?
Here is the full code. C++:
#include<math.h>
#include<iostream>
#include<random>
#include<vector>
#include<algorithm>
#include<fstream>
#include<omp.h>
#include <iomanip>
#include <assert.h>
// parameters
constexpr double v_Rest = -65.0;
constexpr double gNa = 1200.0;
constexpr double gK = 360.0;
constexpr double gL = 3.0;
constexpr double vNa = 115.0;
constexpr double vK = -12.0;
constexpr double vL = 10.6;
constexpr double c = 1.0;
constexpr double knoise = 0.0005;
bool print = false;
bool bisection = false;
bool test = true;
// stepsize PFs
constexpr int steps = 5;
double store[steps];
int prob[steps];
double step[steps];
// time constants
constexpr double t_end = 1.0;
constexpr double delay = 0.1;
constexpr double duration = 0.1;
constexpr double dt = 0.0025;
constexpr int t_steps = t_end/dt;
constexpr int runs = 1000;
double voltage[t_steps];
double doubleRand() {
thread_local std::mt19937 engine(std::random_device{}());
std::normal_distribution<double> distribution(0.0, 1.0);
return distribution(engine);
}
double alphaM(const double v){ return 12.0 * ((2.5 - 0.1 * (v)) / (exp(2.5 - 0.1 * (v)) - 1.0)); }
double betaM(const double v){ return 12.0 * (4.0 * exp(-(v) / 18.0)); }
double betaH(const double v){ return 12.0 * (1.0 / (exp(3.0 - 0.1 * (v)) + 1.0)); }
double alphaH(const double v){ return 12.0 * (0.07 * exp(-(v) / 20.0)); }
double alphaN(const double v){ return 12.0 * ((1.0 - 0.1 * (v)) / (10.0 * (exp(1.0 - 0.1 * (v)) - 1.0))); }
double betaN(const double v){ return 12.0 * (0.125 * exp(-(v) / 80.0)); }
double HH_model(const double I, const double area_factor){
const double A = 1.0e-8 * area_factor;
const double C = c*A;
const double v0 = 0.0;
const double m0 = alphaM(v0)/(alphaM(v0)+betaM(v0));
const double h0 = alphaH(v0)/(alphaH(v0)+betaH(v0));
const double n0 = alphaN(v0)/(alphaN(v0)+betaN(v0));
int count = 0;
for(int j=0; j<runs; j++){
double vT = v0;
double mT = m0;
double hT = h0;
double nT = n0;
for(int i=0; i<t_steps; i++){
double IStim = 0.0;
if ((delay / dt <= (double)i) && ((double)i <= (delay + duration) / dt))
IStim = I;
mT = (mT + dt * alphaM(vT)) / (1.0 + dt * (alphaM(vT) + betaM(vT)));
hT = (hT + dt * alphaH(vT)) / (1.0 + dt * (alphaH(vT) + betaH(vT)));
nT = (nT + dt * alphaN(vT)) / (1.0 + dt * (alphaN(vT) + betaN(vT)));
const double iNa = gNa * pow(mT, 3.0) * hT * (vT - vNa);
const double iK = gK * pow(nT, 4.0) * (vT - vK);
const double iL = gL * (vT-vL);
const double Inoise = (doubleRand() * knoise * sqrt(gNa * A));
const double IIon = ((iNa + iK + iL) * A) + Inoise;
vT += ((-IIon + IStim) / C) * dt;
voltage[i] = vT;
if(vT > 60.0) {
count++;
break;
}
}
}
return count;
}
int main(){
std::cout << HH_model(1.0e-6,1) << std::endl;
}
}
Python:
import matplotlib.pyplot as py
import numpy as np
import scipy.optimize as optimize
from tqdm import tqdm
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
# HH parameters
v_Rest = -65 # in mV
gNa = 1200 # in mS/cm^2
gK = 360 # in mS/cm^2
gL = 0.3*10 # in mS/cm^2
vNa = 115 # in mV
vK = -12 # in mV
vL = 10.6 # in mV
#Number of runs
runs = 1000
c = 1 # in uF/cm^2
def alphaM(v): return 12 * ((2.5 - 0.1 * (v)) / (np.exp(2.5 - 0.1 * (v)) - 1))
def betaM(v): return 12 * (4 * np.exp(-(v) / 18))
def betaH(v): return 12 * (1 / (np.exp(3 - 0.1 * (v)) + 1))
def alphaH(v): return 12 * (0.07 * np.exp(-(v) / 20))
def alphaN(v): return 12 * ((1 - 0.1 * (v)) / (10 * (np.exp(1 - 0.1 * (v)) - 1)))
def betaN(v): return 12 * (0.125 * np.exp(-(v) / 80))
def HH_model(I,area_factor):
count = 0
t_end = 1 # in ms
delay = 0.1 # in ms
duration = 0.1 # in ms
dt = 0.0025 # in ms
area_factor = area_factor
I = I
C = c*A # uF
for j in tqdm(range(0, runs), total=runs):
# Introduction of equations and channels
# compute the timesteps
t_steps= t_end/dt+1
# Compute the initial values
v0 = 0
m0 = alphaM(v0)/(alphaM(v0)+betaM(v0))
h0 = alphaH(v0)/(alphaH(v0)+betaH(v0))
n0 = alphaN(v0)/(alphaN(v0)+betaN(v0))
# Allocate memory for v, m, h, n
v = np.zeros((int(t_steps), 1))
m = np.zeros((int(t_steps), 1))
h = np.zeros((int(t_steps), 1))
n = np.zeros((int(t_steps), 1))
# Set Initial values
v[:, 0] = v0
m[:, 0] = m0
h[:, 0] = h0
n[:, 0] = n0
### Noise component
knoise= 0.0005 #uA/(mS)^1/2
### --------- Step3: SOLVE
for i in range(0, int(t_steps)-1, 1):
# Get current states
vT = v[i]
mT = m[i]
hT = h[i]
nT = n[i]
# Stimulus current
IStim = 0
if delay / dt <= i <= (delay + duration) / dt:
IStim = I # in uA
else:
IStim = 0
# Compute change of m, h and n
m[i + 1] = (mT + dt * alphaM(vT)) / (1 + dt * (alphaM(vT) + betaM(vT)))
h[i + 1] = (hT + dt * alphaH(vT)) / (1 + dt * (alphaH(vT) + betaH(vT)))
n[i + 1] = (nT + dt * alphaN(vT)) / (1 + dt * (alphaN(vT) + betaN(vT)))
# Ionic currents
iNa = gNa * m[i + 1] ** 3. * h[i + 1] * (vT - vNa)
iK = gK * n[i + 1] ** 4. * (vT - vK)
iL = gL * (vT-vL)
Inoise = (np.random.normal(0, 1) * knoise * np.sqrt(gNa * A))
IIon = ((iNa + iK + iL) * A) + Inoise #
# Compute change of voltage
v[i + 1] = (vT + ((-IIon + IStim) / C) * dt)[0] # in ((uA / cm ^ 2) / (uF / cm ^ 2)) * ms == mV
# adjust the voltage to the resting potential
v = v + v_Rest
# test if there was a spike
if max(v[:]-v_Rest) > 60:
count += 1
return count
from Random number generators C++ and Python
No comments:
Post a Comment