计算神经学程序，这是真实神经网络的

#include <stdio.h> #include <stdlib.h> #include <math.h> #include <iostream> //---------------CONSTANTs initialization---------------------------- //--------------- Network Geometry ------------------------------------ #define I_TC 1 // 0 - No layer, 1 - Add Layer #define I_RE 1 // 0 - No layer, 1 - Add Layer #define I_CX 1 // 0 - No layer, 1 - Add Layer #define I_IN 1 // 0 - No layer, 1 - Add Layer #define I_GB 0 // 0 - No GABAb from IN to CX, 1 - Yes #define M 16//32//50 #define Mcx 32//64//100 #define Min 16//32//25 #define M1 16//32 #define Mcx1 32//64 #define Min1 16//32 #define Max ((M > Mcx) ? M : Mcx) #define Max1 ((M1 > Mcx1) ? M1 : Mcx1) //-------------Boundary conditions --------------------------------------- #define BOUND 0 #define BOUNDcx 0 #define SELFING 0 #define SELFINGcx 0 //-------------- Define the connections between cells ----------------- #define MS_RE_RE 5 #define MS_RE_RE1 5 #define MS_RE_RE_MAX ((MS_RE_RE > MS_RE_RE1) ? MS_RE_RE : MS_RE_RE1) #define N_RE_RE (2*MS_RE_RE+1)*(2*MS_RE_RE1+1) #define MS_RE_TC 5 #define MS_RE_TC1 5 #define MS_RE_TC_MAX ((MS_RE_TC > MS_RE_TC1) ? MS_RE_TC : MS_RE_TC1) #define N_RE_TC (2*MS_RE_TC+1)*(2*MS_RE_TC1+1) #define MS_TC_RE 5 #define MS_TC_RE1 5 #define MS_TC_RE_MAX ((MS_TC_RE > MS_TC_RE1) ? MS_TC_RE : MS_TC_RE1) #define N_TC_RE (2*MS_TC_RE+1)*(2*MS_TC_RE1+1) #define MS_CX_CX 5 #define MS_CX_CX1 5 #define MS_CX_CX_MAX ((MS_CX_CX > MS_CX_CX1) ? MS_CX_CX : MS_CX_CX1) #define N_CX_CX (2*MS_CX_CX+1)*(2*MS_CX_CX1+1) #define MS_CX_CX_NMDA 5 #define MS_CX_CX_NMDA1 5 #define MS_CX_CX_NMDA_MAX ((MS_CX_CX_NMDA > MS_CX_CX_NMDA1) ? MS_CX_CX_NMDA : MS_CX_CX_NMDA1) #define N_CX_CX_NMDA (2*MS_CX_CX_NMDA+1)*(2*MS_CX_CX_NMDA1+1) #define MS_CX_IN_NMDA 1 #define MS_CX_IN_NMDA1 1 #define MS_CX_IN_NMDA_MAX ((MS_CX_IN_NMDA > MS_CX_IN_NMDA1) ? MS_CX_IN_NMDA : MS_CX_IN_NMDA1) #define N_CX_IN_NMDA (2*MS_CX_IN_NMDA+1)*Mcx/Min *(2*MS_CX_IN_NMDA1+1)*Mcx1/Min1 #define MS_CX_IN 1 #define MS_CX_IN1 1 #define MS_CX_IN_MAX ((MS_CX_IN > MS_CX_IN1) ? MS_CX_IN : MS_CX_IN1) #define N_CX_IN (2*MS_CX_IN+1)*Mcx/Min*(2*MS_CX_IN1+1)*Mcx1/Min1 #define MS_IN_CX 5 #define MS_IN_CX1 5 #define MS_IN_CX_MAX ((MS_IN_CX > MS_IN_CX1) ? MS_IN_CX : MS_IN_CX1) #define N_IN_CX ((2*MS_IN_CX+1)+4)*Min/Mcx* (2*MS_IN_CX1+4) *Min1/Mcx1 #define MS_TC_CX 10 #define MS_TC_CX1 10 #define MS_TC_CX_MAX ((MS_TC_CX > MS_TC_CX1) ? MS_TC_CX : MS_TC_CX1) #define N_TC_CX (2*MS_TC_CX+1)*M/Mcx*(2*MS_TC_CX1+1) * M1/Mcx1 #define MS_TC_IN 2 #define MS_TC_IN1 2 #define MS_TC_IN_MAX ((MS_TC_IN > MS_TC_IN1) ? MS_TC_IN : MS_TC_IN1) #define N_TC_IN (2*MS_TC_IN+1)* M/Min*(2*MS_TC_IN1+1) * M1/Min1 #define MS_CX_TC 5 #define MS_CX_TC1 5 #define MS_CX_TC_MAX ((MS_CX_TC > MS_CX_TC1) ? MS_CX_TC : MS_CX_TC1) #define N_CX_TC (2*MS_CX_TC+1)* Mcx/M*(2*MS_CX_TC1+1) * Mcx1/M1 //------------Number of ODE for each cell ------------------------------- #define N_RE 7 #define N_TC 12 #define N_GB 2 #define N_DEND 9 #define N_SOMA 4 //3 #define N_CX (N_DEND + N_SOMA) #define N_IN N_CX //4 #define N_EQ1 (N_RE*I_RE + N_TC*I_TC + N_RE_TC*N_GB*I_RE*I_TC)*M*M1 #define N_EQ2 (N_CX*I_CX + N_IN_CX*N_GB*I_CX*I_IN*I_GB)*Mcx*Mcx1 #define N_EQ3 (N_IN*I_IN)*Min*Min1 #define N_EQ (N_EQ1 + N_EQ2 + N_EQ3) // Complete number of ODE //++++++++++++++ Current DESCRIPTION ++++++++++++++++++++++++++++++++++++++ //---------------Low-threshold Ca2+ current (RE cell)--------------------- class IT_RE { static double Shift, Ca_0, Cels; double m_inf, tau_m, h_inf, tau_h, ratio, eca, Phi_m, Phi_h, eca0; public: double iT, m0, h0, Qm, Qh; double G_Ca; IT_RE(double v) { G_Ca = 1.75; Qm = 5; //2.5; //3; Qh = 3; //2.5; //5; Phi_m = pow(Qm,((Cels-24)/10)); Phi_h = pow(Qh,((Cels-24)/10)); m0 = 1/(1 + exp(-(v + Shift + 50)/7.4)); h0 = 1/(1 + exp((v + Shift + 78)/5)); eca0 = 1000*8.31441*(273.15 + Cels)/(2*96489); } void calc(double m, double h, double &fm, double &fh, double v, double cai, double x); }; double IT_RE::Shift = 2, IT_RE::Ca_0 = 2, IT_RE::Cels = 36; void IT_RE::calc(double m, double h, double &fm, double &fh, double v, double cai, double x) { ratio = Ca_0/cai; if(ratio <= 0.)printf("\n LOG ERROR: RE: cai=%lf ratio=%lf",cai,ratio); eca = eca0 * log(ratio); iT = G_Ca*m*m*h*(v - eca); m_inf = 1/(1 + exp(-(v + 52)/7.4)); tau_m = (3 + 1/(exp((v + 27)/10) + exp(-(v + 102)/15)))/Phi_m; h_inf = 1/(1 + exp((v + 80)/5)); tau_h = (85 + 1/(exp((v + 48)/4) + exp(-(v + 407)/50)))/Phi_h; fm = -(m - m_inf)/tau_m; fh = -(h - h_inf)/tau_h; } //--------------fast Na and K current (RE and TC cells)------------------ class INaK { static double Cels; double Alpha1, Beta1, Alpha2, Beta2, Alpha3, Beta3, v2, v2K, Phi; double tau_m, m_inf, tau_h, h_inf, tau_n, n_inf; public: static double E_Na, E_K; double iK, iNa, m0, h0, n0; double G_Na, G_K, Vtr, VtrK; INaK(double v) { G_K = 10;/////////////////////// G_Na = 100;///////////////////// Vtr = -50; VtrK = -50; v2 = v - Vtr; v2K = v - VtrK; Phi = pow(3,((Cels-36)/10)); Alpha1 = 0.32*(13 - v2)/(exp((13 - v2)/4) - 1); Beta1 = 0.28*(v2 - 40)/(exp((v2 - 40)/5) - 1); m0 = Alpha1/(Alpha1 + Beta1); Alpha2 = 0.128*exp((17 - v2)/18); Beta2 = 4/(exp((40 - v2)/5) + 1); h0 = Alpha2/(Alpha2 + Beta2); Alpha3 = 0.032*(15 - v2K)/(exp((15 - v2K)/5) - 1); Beta3 = 0.5*exp((10 - v2K)/40); n0 = Alpha3/(Alpha3 + Beta3); } void calc(double m, double h, double n, double &fm, double &fh, double &fn, double v, double x); }; double INaK::E_K = -95, INaK::E_Na = 50, INaK::Cels = 36; void INaK::calc(double m, double h, double n, double &fm, double &fh, double &fn, double v, double x){ v2 = v - Vtr; v2K = v - VtrK; iNa = G_Na*m*m*m*h*(v - E_Na); Alpha1 = 0.32*(13 - v2)/(exp((13 - v2)/4) - 1); Beta1 = 0.28*(v2 - 40)/(exp((v2 - 40)/5) - 1); tau_m = 1/(Alpha1 + Beta1) / Phi; m_inf = Alpha1/(Alpha1 + Beta1); Alpha2 = 0.128*exp((17 - v2)/18); Beta2 = 4/(exp((40 - v2)/5) + 1); tau_h = 1/(Alpha2 + Beta2) / Phi; h_inf = Alpha2/(Alpha2 + Beta2); fm = -(m - m_inf)/tau_m; fh = -(h - h_inf)/tau_h; iK = G_K* n*n*n*n*(v - E_K); Alpha3 = 0.032*(15 - v2K)/(exp((15 - v2K)/5) - 1); Beta3 = 0.5*exp((10 - v2K)/40); tau_n = 1/(Alpha3 + Beta3) / Phi; n_inf = Alpha3/(Alpha3 + Beta3); fn = -(n - n_inf)/tau_n; } //------------------Ca-dynamics------------------------------------ class ICa { static double Ca_inf, K_T, K_d; double drive, drive0; public: double Taur, D; ICa() {Taur = 5; D = 1.; //0.1; drive0 = 10.0/(2.*96489.); } void calc(double cai, double &fcai, double iT, double x); }; double ICa::Ca_inf = 2.4e-4; double ICa::K_T = 0.0001, ICa::K_d = 0.0001; void ICa::calc(double cai, double &fcai, double iT, double x) { drive = -drive0 * iT / D; if(drive < 0) drive = 0; fcai = drive + (Ca_inf - cai)/Taur; // - K_T*cai/(cai + K_d); } //------------------Low-theshold Ca2+ current (TC cell)----------------- class IT_TC { static double Ca_0, Cels, Qm, Qh, Shift; double m_inf, tau_m, h_inf, tau_h, Phi_h, Phi_m, ratio, eca, eca0; public: double iT, m0, h0; double G_Ca; IT_TC(double v) { G_Ca = 2; //2;/////////////////////////////////// Phi_m = pow(Qm,((Cels-24)/10)); Phi_h = pow(Qh,((Cels-24)/10)); m0 = 1 / (1+exp(-(v+59)/6.2));//////////////////////// h0 = 1 / (1+exp