我已经用C++和Python编程了一个模型。这个模型有一个嘈杂的输入组件，我可以用这个C++：来代替它
double doubleRand() {
thread_local std::mt19937 generator(std::random_device{}());
std::normal_distribution<double> distribution(0.0, 1.0);
return distribution(generator);
}
或者这个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]
以下内容非常奇怪：如果我省略了有噪声的组件(Inoise=0(，那么两个模型(C++和Python(都会给出完全相同的结果。如果我只引入噪声分量(Istim=0(，那么两个模型都会给出结果(即，在1000次运行时几乎没有差异的自然波动(。但是，如果我选择Istim=0.000001并添加噪波，则结果会相差30%。这怎么可能？
这是完整的代码。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
您在Python代码中搞砸了缩进。这些线路
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)))
条件delay / dt <= i <= (delay + duration) / dt为True 时不执行
在修复缩进后，Python代码生成866，这几乎与C++代码的876结果相匹配。