/***************************************************************************
* LIFTimeDrivenModel_1_2.cpp *
* ------------------- *
* copyright : (C) 2013 by Jesus Garrido and Francisco Naveros *
* email : jgarrido@atc.ugr.es, 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/TimeDrivenPurkinjeCell.h"
#include "../../include/neuron_model/VectorNeuronState.h"
#include <iostream>
#include <cmath>
#include <string>
#include "../../include/openmp/openmp.h"
#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"
#include "../../include/openmp/openmp.h"
void TimeDrivenPurkinjeCell::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->eexc)==1){
// skip_comments(fh,Currentline);
// if (fscanf(fh,"%f",&this->einh)==1){
// skip_comments(fh,Currentline);
// if(fscanf(fh,"%f",&this->erest)==1){
// skip_comments(fh,Currentline);
// if(fscanf(fh,"%f",&this->vthr)==1){
// skip_comments(fh,Currentline);
// if(fscanf(fh,"%f",&this->cm)==1){
// inv_cm=1.0f/cm;
// skip_comments(fh,Currentline);
// if(fscanf(fh,"%f",&this->texc)==1){
// inv_texc=1.0f/texc;
// skip_comments(fh,Currentline);
// if(fscanf(fh,"%f",&this->tinh)==1){
// inv_tinh=1.0f/tinh;
// skip_comments(fh,Currentline);
// if(fscanf(fh,"%f",&this->tref)==1){
// skip_comments(fh,Currentline);
// if(fscanf(fh,"%f",&this->grest)==1){
// skip_comments(fh,Currentline);
this->InitialState = (VectorNeuronState *) new VectorNeuronState(N_NeuronStateVariables, true);
// } 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 TimeDrivenPurkinjeCell::SynapsisEffect(int index, Interconnection * InputConnection){
this->GetVectorNeuronState()->IncrementStateVariableAtCPU(index,N_DifferentialNeuronState+InputConnection->GetType(),1e-6f*InputConnection->GetWeight());
}
TimeDrivenPurkinjeCell::TimeDrivenPurkinjeCell(string NeuronTypeID, string NeuronModelID): TimeDrivenNeuronModel(NeuronTypeID, NeuronModelID), g_L(0.02f),
g_Ca(0.001f), g_M(0.75f), Cylinder_length_of_the_soma(0.0015f), Radius_of_the_soma(0.0008f), Area(3.141592f*0.0015f*2.0f*0.0008f),
inv_Area(1.0f/(3.141592f*0.0015f*2.0f*0.0008f)), Membrane_capacitance(1.0f), inv_Membrane_capacitance(1.0f/1.0f)
{
eexc=0.0f;
einh=-80.0f ;
vthr= -35.0f;
erest=-65.0f;
texc=1.0f;
inv_texc=1.0f/texc;
tinh=2;
inv_tinh=1.0f/tinh;
tref=1.35f;
tref_0_5=tref*0.5f;
inv_tref_0_5=1.0f/tref_0_5;
spkpeak=31.0f;
}
TimeDrivenPurkinjeCell::~TimeDrivenPurkinjeCell(void)
{
}
void TimeDrivenPurkinjeCell::LoadNeuronModel() throw (EDLUTFileException){
this->LoadNeuronModel(this->GetModelID()+".cfg");
}
VectorNeuronState * TimeDrivenPurkinjeCell::InitializeState(){
return this->GetVectorNeuronState();
}
InternalSpike * TimeDrivenPurkinjeCell::ProcessInputSpike(Interconnection * inter, Neuron * target, double time){
// Add the effect of the input spike
this->SynapsisEffect(target->GetIndex_VectorNeuronState(),inter);
return 0;
}
bool TimeDrivenPurkinjeCell::UpdateState(int index, VectorNeuronState * State, double CurrentTime){
//float * NeuronState;
//NeuronState[0] --> V
//NeuronState[1] --> c
//NeuronState[2] --> M
//NeuronState[3] --> gexc
//NeuronState[4] --> ginh
bool * internalSpike=State->getInternalSpike();
double last_update = State->GetLastUpdateTime(0);
double elapsed_time = CurrentTime - last_update;
float elapsed_time_f=elapsed_time;
for(int j=0; j<NumberOfOpenMPTasks-1; j++){
#ifdef _OPENMP
#if _OPENMP >= OPENMPVERSION30
#pragma omp task firstprivate (j) shared(internalSpike, State, CurrentTime)
#endif
#endif
{
for (int i=LimitOfOpenMPTasks[j]; i< LimitOfOpenMPTasks[j+1]; i++){
State->AddElapsedTime(i,elapsed_time);
double last_spike = State->GetLastSpikeTime(i);
last_spike*=1000;//ms
float * NeuronState=State->GetStateVariableAt(i);
bool spike = false;
NeuronState[5]=last_spike;
this->integrationMethod->NextDifferentialEcuationValue(i, NeuronState, elapsed_time_f*1000);
if (last_spike > this->tref && NeuronState[0] > this->vthr) {
State->NewFiredSpike(i);
last_spike=0;
spike = true;
}
if(last_spike < tref){
if(last_spike <= tref_0_5){
NeuronState[0]=vthr+(spkpeak-vthr)*(last_spike*inv_tref_0_5);
}else{
NeuronState[0]=spkpeak-(spkpeak-erest)*((last_spike-tref_0_5)*inv_tref_0_5);
}
}
internalSpike[i]=spike;
State->SetLastUpdateTime(i,CurrentTime);
}
}
}
for (int i=LimitOfOpenMPTasks[NumberOfOpenMPTasks-1]; i< LimitOfOpenMPTasks[NumberOfOpenMPTasks]; i++){
State->AddElapsedTime(i,elapsed_time);
double last_spike = State->GetLastSpikeTime(i);
last_spike*=1000;//ms
float * NeuronState=State->GetStateVariableAt(i);
bool spike = false;
NeuronState[5]=last_spike;
this->integrationMethod->NextDifferentialEcuationValue(i, NeuronState, elapsed_time_f*1000);
if (last_spike > this->tref && NeuronState[0] > this->vthr) {
State->NewFiredSpike(i);
last_spike=0;
spike = true;
}
if(last_spike < tref){
if(last_spike <= tref_0_5){
NeuronState[0]=vthr+(spkpeak-vthr)*(last_spike*inv_tref_0_5);
}else{
NeuronState[0]=spkpeak-(spkpeak-erest)*((last_spike-tref_0_5)*inv_tref_0_5);
}
}
internalSpike[i]=spike;
State->SetLastUpdateTime(i,CurrentTime);
}
#ifdef _OPENMP
#if _OPENMP >= OPENMPVERSION30
#pragma omp taskwait
#endif
#endif
return false;
}
ostream & TimeDrivenPurkinjeCell::PrintInfo(ostream & out){
//out << "- Leaky Time-Driven Model: " << this->GetModelID() << endl;
//out << "\tExc. Reversal Potential: " << this->eexc << "V\tInh. Reversal Potential: " << this->einh << "V\tResting potential: " << this->erest << "V" << endl;
//out << "\tFiring threshold: " << this->vthr << "V\tMembrane capacitance: " << this->cm << "nS\tExcitatory Time Constant: " << this->texc << "s" << endl;
//out << "\tInhibitory time constant: " << this->tinh << "s\tRefractory Period: " << this->tref << "s\tResting Conductance: " << this->grest << "nS" << endl;
return out;
}
void TimeDrivenPurkinjeCell::InitializeStates(int N_neurons, int OpenMPQueueIndex){
//Initialize neural state variables.
float alpha_ca, beta_ca, alpha_M, beta_M;
Get_alpha_and_beta_values(erest, &alpha_ca, &beta_ca, &alpha_M, &beta_M);
//c_inf
float c_inf=alpha_ca/(alpha_ca+beta_ca);
//M_inf
float M_inf=alpha_M/(alpha_M+beta_M);
float initialization[] = {erest,c_inf,M_inf,0.0f,0.0f};
InitialState->InitializeStates(N_neurons, initialization);
//Initialize integration method state variables.
this->integrationMethod->InitializeStates(N_neurons, initialization);
//Calculate number of OpenMP task and size of each one.
CalculateTaskSizes(N_neurons, 50);
}
void TimeDrivenPurkinjeCell::EvaluateDifferentialEcuation(float * NeuronState, float * AuxNeuronState){
float V=NeuronState[0];
float ca=NeuronState[1];
float M=NeuronState[2];
float g_exc=NeuronState[3];
float g_inh=NeuronState[4];
float last_spike=NeuronState[5];
//V
if(last_spike >= tref){
AuxNeuronState[0]=(-g_L*(V+70.0f)-g_Ca*ca*ca*(V-125.0f)-g_M*M*(V+95.0f) + (g_exc * (this->eexc - V) + g_inh * (this->einh - V))*inv_Area )*inv_Membrane_capacitance;
}else if(last_spike <= tref_0_5){
AuxNeuronState[0]=(spkpeak-vthr)*inv_tref_0_5;
}else{
AuxNeuronState[0]=(erest-spkpeak)*inv_tref_0_5;
}
float alpha_ca, beta_ca, alpha_M, beta_M;
Get_alpha_and_beta_values(V, &alpha_ca, &beta_ca, &alpha_M, &beta_M);
//ca
float inv_tau_ca=(alpha_ca+beta_ca);
AuxNeuronState[1]=alpha_ca - ca*inv_tau_ca;
//M
float inv_tau_M=(alpha_M+beta_M);
AuxNeuronState[2]=alpha_M - M*inv_tau_M;
}
void TimeDrivenPurkinjeCell::EvaluateTimeDependentEcuation(float * NeuronState, float elapsed_time){
float limit=1e-20;
if(NeuronState[N_DifferentialNeuronState]<limit){
NeuronState[N_DifferentialNeuronState]=0.0f;
}else{
NeuronState[N_DifferentialNeuronState]*= ExponentialTable::GetResult(-(elapsed_time*this->inv_texc));
}
if(NeuronState[N_DifferentialNeuronState+1]<limit){
NeuronState[N_DifferentialNeuronState+1]=0.0f;
}else{
NeuronState[N_DifferentialNeuronState+1]*= ExponentialTable::GetResult(-(elapsed_time*this->inv_tinh));
}
}
int TimeDrivenPurkinjeCell::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;
}
}
const float TimeDrivenPurkinjeCell::Max_V=35.0f;
const float TimeDrivenPurkinjeCell::Min_V=-100.0f;
const float TimeDrivenPurkinjeCell::aux=(TimeDrivenPurkinjeCell::TableSize-1)/( TimeDrivenPurkinjeCell::Max_V- TimeDrivenPurkinjeCell::Min_V);
float * TimeDrivenPurkinjeCell::alpha_ca_table=Generate_alpha_ca_table();
float * TimeDrivenPurkinjeCell::beta_ca_table=Generate_beta_ca_table();
float * TimeDrivenPurkinjeCell::alpha_M_table=Generate_alpha_M_table();
float * TimeDrivenPurkinjeCell::beta_M_table=Generate_beta_M_table();