/***************************************************************************
* EdidioGranuleCell_TimeDriven.cpp *
* ------------------- *
* copyright : (C) 2013 by Francisco Naveros *
* email : fnaveros@atc.ugr.es *
***************************************************************************/
/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 3 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
#include "../../include/neuron_model/EgidioGranuleCell_TimeDriven.h"
#include "../../include/neuron_model/VectorNeuronState.h"
#include <iostream>
#include <cmath>
#include <string>
#include "../../include/openmp/openmp.h"
//This neuron model is implemented in milisecond. EDLUT is implemented in second and it is necesary to
//use this constant in order to adapt this model to EDLUT.
#define ms_to_s 1000.0f
#include "../../include/spike/EDLUTFileException.h"
#include "../../include/spike/Neuron.h"
#include "../../include/spike/InternalSpike.h"
#include "../../include/spike/PropagatedSpike.h"
#include "../../include/spike/Interconnection.h"
#include "../../include/simulation/Utils.h"
void EgidioGranuleCell_TimeDriven::LoadNeuronModel(string ConfigFile) throw (EDLUTFileException){
FILE *fh;
long Currentline = 0L;
fh=fopen(ConfigFile.c_str(),"rt");
if(fh){
Currentline=1L;
skip_comments(fh,Currentline);
if(fscanf(fh,"%f",&this->gMAXNa_f)==1){
skip_comments(fh,Currentline);
if (fscanf(fh,"%f",&this->gMAXNa_r)==1){
skip_comments(fh,Currentline);
if(fscanf(fh,"%f",&this->gMAXNa_p)==1){
skip_comments(fh,Currentline);
if(fscanf(fh,"%f",&this->gMAXK_V)==1){
skip_comments(fh,Currentline);
if(fscanf(fh,"%f",&this->gMAXK_A)==1){
skip_comments(fh,Currentline);
if(fscanf(fh,"%f",&this->gMAXK_IR)==1){
skip_comments(fh,Currentline);
if(fscanf(fh,"%f",&this->gMAXK_Ca)==1){
skip_comments(fh,Currentline);
if(fscanf(fh,"%f",&this->gMAXCa)==1){
skip_comments(fh,Currentline);
if(fscanf(fh,"%f",&this->gMAXK_sl)==1){
skip_comments(fh,Currentline);
this->InitialState = (VectorNeuronState *) new VectorNeuronState(17, true);
}
//TODOS LOS CODIGOS DE ERROR HAY QUE MODIFICARLOS, PORQUE AHORA LOS PARAMETROS QUE SE CARGAN SON OTROS.
else {
throw EDLUTFileException(13,60,3,1,Currentline);
}
} else {
throw EDLUTFileException(13,61,3,1,Currentline);
}
} else {
throw EDLUTFileException(13,62,3,1,Currentline);
}
} else {
throw EDLUTFileException(13,63,3,1,Currentline);
}
} else {
throw EDLUTFileException(13,64,3,1,Currentline);
}
} else {
throw EDLUTFileException(13,65,3,1,Currentline);
}
} else {
throw EDLUTFileException(13,66,3,1,Currentline);
}
} else {
throw EDLUTFileException(13,67,3,1,Currentline);
}
} else {
throw EDLUTFileException(13,68,3,1,Currentline);
}
//INTEGRATION METHOD
this->integrationMethod = LoadIntegrationMethod::loadIntegrationMethod((TimeDrivenNeuronModel *)this, fh, &Currentline, N_NeuronStateVariables, N_DifferentialNeuronState, N_TimeDependentNeuronState);
}
}
void EgidioGranuleCell_TimeDriven::SynapsisEffect(int index, Interconnection * InputConnection){
this->GetVectorNeuronState()->IncrementStateVariableAtCPU(index,N_DifferentialNeuronState+InputConnection->GetType(),1e-9f*InputConnection->GetWeight());
}
EgidioGranuleCell_TimeDriven::EgidioGranuleCell_TimeDriven(string NeuronTypeID, string NeuronModelID): TimeDrivenNeuronModel(NeuronTypeID, NeuronModelID), gMAXNa_f(0), gMAXNa_r(0), gMAXNa_p(0), gMAXK_V(0), gMAXK_A(0), gMAXK_IR(0), gMAXK_Ca(0),
gMAXCa(0), gMAXK_sl(0),
gLkg1(5.68e-5),
gLkg2(2.17e-5),
VNa(87.39),
VK(-84.69),
VLkg1(-58),
VLkg2(-65),
V0_xK_Ai(-46.7),
K_xK_Ai(-19.8),
V0_yK_Ai(-78.8),
K_yK_Ai(8.4),
V0_xK_sli(-30),
B_xK_sli(6),
F(96485.309),
A(1e-04),
d(0.2),
betaCa(1.5),
Ca0(1e-04),
R(8.3134),
cao(2),
Cm(1.0e-3),
temper(30),
Q10_20 ( pow(3,((temper-20)/10))),
Q10_22 ( pow(3,((temper-22)/10))),
Q10_30 ( pow(3,((temper-30)/10))),
Q10_6_3 ( pow(3,((temper-6.3)/10))),
//This is a constant current which can be externally injected to the cell.
I_inj_abs(11e-12)/*I_inj_abs(0)*/,
I_inj(-I_inj_abs*1000/299.26058e-8),
eexc(0.0),
einh(-80),
texc(0.5),
tinh(10),
vthr(-0.25)
{
}
EgidioGranuleCell_TimeDriven::~EgidioGranuleCell_TimeDriven(void)
{
}
void EgidioGranuleCell_TimeDriven::LoadNeuronModel() throw (EDLUTFileException){
this->LoadNeuronModel(this->GetModelID()+".cfg");
}
VectorNeuronState * EgidioGranuleCell_TimeDriven::InitializeState(){
return this->GetVectorNeuronState();
}
InternalSpike * EgidioGranuleCell_TimeDriven::ProcessInputSpike(Interconnection * inter, Neuron * target, double time){
// Add the effect of the input spike
this->SynapsisEffect(target->GetIndex_VectorNeuronState(),inter);
return 0;
}
float EgidioGranuleCell_TimeDriven::nernst(float ci, float co, float z, float temper){
//return (1000*(R*(temper + 273.15f)/F)/z*log(co/ci));
return (1000*(R*(temper + 273.15f)/F)/z*log(abs(co/ci)));
}
float EgidioGranuleCell_TimeDriven::linoid(float x, float y){
float f=0.0;
if (abs(x/y)<1e-06f){
f=y*(1-x/y/2);
}else{
f=x/(ExponentialTable::GetResult(x/y)-1);
}
return f;
}
bool EgidioGranuleCell_TimeDriven::UpdateState(int index, VectorNeuronState * State, double CurrentTime){
bool * internalSpike=State->getInternalSpike();
int Size=State->GetSizeState();
double last_update;
double elapsed_time;
float elapsed_time_f;
double last_spike;
bool spike;
float vm_cou;
int i;
float previous_V;
float * NeuronState;
for (int i=0; i< Size; i++){
last_update = State->GetLastUpdateTime(i);
elapsed_time = CurrentTime - last_update;
elapsed_time_f=elapsed_time;
State->AddElapsedTime(i,elapsed_time);
last_spike = State->GetLastSpikeTime(i);
NeuronState=State->GetStateVariableAt(i);
spike = false;
previous_V=NeuronState[14];
this->integrationMethod->NextDifferentialEcuationValue(i, NeuronState, elapsed_time_f);
if(NeuronState[14]>vthr && previous_V<vthr){
State->NewFiredSpike(i);
spike = true;
}
internalSpike[i]=spike;
State->SetLastUpdateTime(i,CurrentTime);
}
return false;
}
ostream & EgidioGranuleCell_TimeDriven::PrintInfo(ostream & out){
return out;
}
void EgidioGranuleCell_TimeDriven::InitializeStates(int N_neurons, int OpenMPQueueIndex){
//Initial State
float xNa_f=0.00047309535f;
float yNa_f=1.0f;
float xNa_r=0.00013423511f;
float yNa_r=0.96227829f;
float xNa_p=0.00050020111f;
float xK_V=0.010183001f;
float xK_A=0.15685486f;
float yK_A=0.53565367f;
float xK_IR=0.37337035f;
float xK_Ca=0.00012384122f;
float xCa=0.0021951104f;
float yCa=0.89509747f;
float xK_sl=0.00024031171f;
float Ca=Ca0;
float V=-80.0f;
float gexc=0.0f;
float ginh=0.0f;
//Initialize neural state variables.
float initialization[] = {xNa_f,yNa_f,xNa_r,yNa_r,xNa_p,xK_V,xK_A,yK_A,xK_IR,xK_Ca,xCa,yCa,xK_sl,Ca,V,gexc,ginh};
InitialState->InitializeStates(N_neurons, initialization);
//Initialize integration method state variables.
this->integrationMethod->InitializeStates(N_neurons, initialization);
}
void EgidioGranuleCell_TimeDriven::EvaluateDifferentialEcuation(float * NeuronState, float * AuxNeuronState){
float previous_V=NeuronState[14];
float VCa=nernst(NeuronState[13],cao,2,temper);
float alphaxNa_f = Q10_20*(-0.3f)*linoid(previous_V+19, -10);
float betaxNa_f = Q10_20*12*ExponentialTable::GetResult(-(previous_V+44)/18.182f);
float xNa_f_inf = alphaxNa_f/(alphaxNa_f + betaxNa_f);
float inv_tauxNa_f = (alphaxNa_f + betaxNa_f);
float alphayNa_f = Q10_20*0.105f*ExponentialTable::GetResult(-(previous_V+44)/3.333f);
float betayNa_f = Q10_20*1.5f/(1+ExponentialTable::GetResult(-(previous_V+11)/5));
float yNa_f_inf = alphayNa_f/(alphayNa_f + betayNa_f);
float inv_tauyNa_f = (alphayNa_f + betayNa_f);
float alphaxNa_r = Q10_20*(0.00008f-0.00493f*linoid(previous_V-4.48754f,-6.81881f));
float betaxNa_r = Q10_20*(0.04752f+0.01558f*linoid(previous_V+43.97494f,0.10818f));
float xNa_r_inf = alphaxNa_r/(alphaxNa_r + betaxNa_r);
float inv_tauxNa_r = (alphaxNa_r + betaxNa_r);
float alphayNa_r = Q10_20*0.31836f*ExponentialTable::GetResult(-(previous_V+80)/62.52621f);
float betayNa_r = Q10_20*0.01014f*ExponentialTable::GetResult((previous_V+83.3332f)/16.05379f);
float yNa_r_inf = alphayNa_r/(alphayNa_r + betayNa_r);
float inv_tauyNa_r = (alphayNa_r + betayNa_r);
float alphaxNa_p = Q10_30*(-0.091f)*linoid(previous_V+42,-5);
float betaxNa_p = Q10_30*0.062f*linoid(previous_V+42,5);
float xNa_p_inf = 1/(1+ExponentialTable::GetResult(-(previous_V+42)/5));
float inv_tauxNa_p = (alphaxNa_p + betaxNa_p)*0.2f;
float alphaxK_V = Q10_6_3*(-0.01f)*linoid(previous_V+25,-10);
float betaxK_V = Q10_6_3*0.125f*ExponentialTable::GetResult(-0.0125f*(previous_V+35));
float xK_V_inf = alphaxK_V/(alphaxK_V + betaxK_V);
float inv_tauxK_V = (alphaxK_V + betaxK_V);
float alphaxK_A = (Q10_20*4.88826f)/(1+ExponentialTable::GetResult(-(previous_V+9.17203f)/23.32708f));
float betaxK_A = (Q10_20*0.99285f)/ExponentialTable::GetResult((previous_V+18.27914f)/19.47175f);
float xK_A_inf = 1/(1+ExponentialTable::GetResult((previous_V-V0_xK_Ai)/K_xK_Ai));
float inv_tauxK_A = (alphaxK_A + betaxK_A);
float alphayK_A = (Q10_20*0.11042f)/(1+ExponentialTable::GetResult((previous_V+111.33209f)/12.8433f));
float betayK_A = (Q10_20*0.10353f)/(1+ExponentialTable::GetResult(-(previous_V+49.9537f)/8.90123f));
float yK_A_inf = 1/(1+ExponentialTable::GetResult((previous_V-V0_yK_Ai)/K_yK_Ai));
float inv_tauyK_A = (alphayK_A + betayK_A);
float alphaxK_IR = Q10_20*0.13289f*ExponentialTable::GetResult(-(previous_V+83.94f)/24.3902f);
float betaxK_IR = Q10_20*0.16994f*ExponentialTable::GetResult((previous_V+83.94f)/35.714f);
float xK_IR_inf = alphaxK_IR/(alphaxK_IR + betaxK_IR);
float inv_tauxK_IR = (alphaxK_IR + betaxK_IR);
float alphaxK_Ca = (Q10_30*2.5f)/(1+(0.0015f*ExponentialTable::GetResult(-previous_V/11.765f))/NeuronState[13]);
float betaxK_Ca = (Q10_30*1.5f)/(1+NeuronState[13]/(0.00015*ExponentialTable::GetResult(-previous_V/11.765f)));
float xK_Ca_inf = alphaxK_Ca/(alphaxK_Ca + betaxK_Ca);
float inv_tauxK_Ca = (alphaxK_Ca + betaxK_Ca);
float alphaxCa = Q10_20*0.04944f*ExponentialTable::GetResult((previous_V+29.06f)/15.87301587302f);
float betaxCa = Q10_20*0.08298f*ExponentialTable::GetResult(-(previous_V+18.66f)/25.641f);
float xCa_inf = alphaxCa/(alphaxCa + betaxCa);
float inv_tauxCa = (alphaxCa + betaxCa);
float alphayCa = Q10_20*0.0013f*ExponentialTable::GetResult(-(previous_V+48)/18.183f);
float betayCa = Q10_20*0.0013f*ExponentialTable::GetResult((previous_V+48)/83.33f);
float yCa_inf = alphayCa/(alphayCa + betayCa);
float inv_tauyCa = (alphayCa + betayCa);
float alphaxK_sl = Q10_22*0.0033f*ExponentialTable::GetResult((previous_V+30)/40);
float betaxK_sl = Q10_22*0.0033f*ExponentialTable::GetResult(-(previous_V+30)/20);
float xK_sl_inf = 1/(1+ExponentialTable::GetResult(-(previous_V-V0_xK_sli)/B_xK_sli));
float inv_tauxK_sl = (alphaxK_sl + betaxK_sl);
float gNa_f = gMAXNa_f * NeuronState[0]*NeuronState[0]*NeuronState[0] * NeuronState[1];
float gNa_r = gMAXNa_r * NeuronState[2] * NeuronState[3];
float gNa_p= gMAXNa_p * NeuronState[4];
float gK_V = gMAXK_V * NeuronState[5]*NeuronState[5]*NeuronState[5]*NeuronState[5];
float gK_A = gMAXK_A * NeuronState[6]*NeuronState[6]*NeuronState[6] * NeuronState[7];
float gK_IR = gMAXK_IR * NeuronState[8];
float gK_Ca=gMAXK_Ca * NeuronState[9];
float gCa = gMAXCa * NeuronState[10]*NeuronState[10] * NeuronState[11];
float gK_sl = gMAXK_sl * NeuronState[12];
AuxNeuronState[0]=ms_to_s*(xNa_f_inf - NeuronState[0])*inv_tauxNa_f;
AuxNeuronState[1]=ms_to_s*(yNa_f_inf - NeuronState[1])*inv_tauyNa_f;
AuxNeuronState[2]=ms_to_s*(xNa_r_inf - NeuronState[2])*inv_tauxNa_r;
AuxNeuronState[3]=ms_to_s*(yNa_r_inf - NeuronState[3])*inv_tauyNa_r;
AuxNeuronState[4]=ms_to_s*(xNa_p_inf - NeuronState[4])*inv_tauxNa_p;
AuxNeuronState[5]=ms_to_s*(xK_V_inf - NeuronState[5])*inv_tauxK_V;
AuxNeuronState[6]=ms_to_s*(xK_A_inf - NeuronState[6])*inv_tauxK_A;
AuxNeuronState[7]=ms_to_s*(yK_A_inf - NeuronState[7])*inv_tauyK_A;
AuxNeuronState[8]=ms_to_s*(xK_IR_inf - NeuronState[8])*inv_tauxK_IR;
AuxNeuronState[9]=ms_to_s*(xK_Ca_inf - NeuronState[9])*inv_tauxK_Ca;
AuxNeuronState[10]=ms_to_s*(xCa_inf - NeuronState[10])*inv_tauxCa;
AuxNeuronState[11]=ms_to_s*(yCa_inf - NeuronState[11])*inv_tauyCa;
AuxNeuronState[12]=ms_to_s*(xK_sl_inf-NeuronState[12])*inv_tauxK_sl;
AuxNeuronState[13]=ms_to_s*(-gCa*(previous_V-VCa)/(2*F*A*d) - (betaCa*(NeuronState[13] - Ca0)));
AuxNeuronState[14]=ms_to_s*(-1/Cm)*((NeuronState[15]/299.26058e-8f) * (previous_V - eexc) + (NeuronState[16]/299.26058e-8f) * (previous_V - einh)+gNa_f*(previous_V-VNa)+gNa_r*(previous_V-VNa)+gNa_p*(previous_V-VNa)+gK_V*(previous_V-VK)+gK_A*(previous_V-VK)+gK_IR*(previous_V-VK)+gK_Ca*(previous_V-VK)+gCa*(previous_V-VCa)+gK_sl*(previous_V-VK)+gLkg1*(previous_V-VLkg1)+gLkg2*(previous_V-VLkg2)+I_inj);
}
void EgidioGranuleCell_TimeDriven::EvaluateTimeDependentEcuation(float * NeuronState, float elapsed_time){
//NeuronState[15]*= ExponentialTable::GetResult(-(ms_to_s*elapsed_time/this->texc));
//NeuronState[16]*= ExponentialTable::GetResult(-(ms_to_s*elapsed_time/this->tinh));
if(NeuronState[15]<1e-30){
NeuronState[15]=0.0f;
}else{
NeuronState[15]*= ExponentialTable::GetResult(-(ms_to_s*elapsed_time/this->texc));
}
if(NeuronState[16]<1e-30){
NeuronState[16]=0.0f;
}else{
NeuronState[16]*= ExponentialTable::GetResult(-(ms_to_s*elapsed_time/this->tinh));
}
}
int EgidioGranuleCell_TimeDriven::CheckSynapseTypeNumber(int Type){
if(Type<N_TimeDependentNeuronState && Type>=0){
return Type;
}else{
cout<<"Neuron model "<<this->GetTypeID()<<", "<<this->GetModelID()<<" does not support input synapses of type "<<Type<<endl;
return 0;
}
}