Neurite: electrophysiological-mechanical coupling simulation framework (Garcia-Grajales et al 2015)

 Download zip file 
Help downloading and running models
Accession:168861
Neurite simulates the electrical signal propagation in myelinated and unmyelinated axons, and in dendritic trees under mechanical loading. Two different solvers are available (explicit and implicit) with sequential (CPU) and parallel (GPUs) versions
Reference:
1 . García-Grajales JA, Rucabado G, García-Dopico A, Peña JM, Jérusalem A (2015) Neurite, a finite difference large scale parallel program for the simulation of electrical signal propagation in neurites under mechanical loading. PLoS One 10:e0116532 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Axon; Dendrite;
Brain Region(s)/Organism:
Cell Type(s): Myelinated neuron;
Channel(s): I Sodium; I Potassium;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: C or C++ program;
Model Concept(s): Action Potential Initiation; Axonal Action Potentials; Action Potentials;
Implementer(s): Garcia-Grajales, Julian Andres ;
Search NeuronDB for information about:  I Sodium; I Potassium;
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////// $$    $  $$$$  $   $  $$$$$  $ $$$$$$$ $$$$  ///////////////////////////
/////////////////////////// $ $   $  $     $   $  $   $  $    $    $     ///////////////////////////
/////////////////////////// $  $  $  $$$$  $   $  $$$$$  $    $    $$$$  ///////////////////////////
/////////////////////////// $   $ $  $     $   $  $  $   $    $    $     ///////////////////////////
/////////////////////////// $    $$  $$$$  $$$$$  $   $  $    $    $$$$  ///////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////

/*!
 * \mainpage NEURITE
 * \brief Neurite electrical signal propagation under mechanical loading in sequential and parallel simulations.
 * \author Julián Andrés García Grajales
 * \author Jose María Peña
 * \author Antoine Jérusalem
 * \copyright (C) 2014  This software is licenced under the terms stipulated in the license.txt file located in its root directory
 */

//
//
// File author(s):  <Julian Andres Garcia Grajales, Jose Maria Pena and Antoine Jerusalem>, (C) 2014
//
// Copyright: this software is licenced under the terms stipulated in the license.txt file located in its root directory
//
//

/*!
 * \file Neurite.cpp
 * \brief This is the main file for Neurite.
 * \details The time looping is in this file.
 */
#include "configuration.h"
#include "neurite.h"  
#include "graphics.h"
#include "init.h"
#include "input.h"
#include "coupling.h"
#include "discretization.h"
#include "mechanical_model.h"
#include "space.h"
#include "GL/glut.h"
#include <time.h>
#include "solver.h"

#if defined GPU
#include "cudaException.h"
#include "solverEx-gpu.h"
#include "solverIm-gpu.h"

#elif defined CPU
#include "solverEx-cpu.h"
#include "solverIm-cpu.h"
#else
#error "You did not choose any processor"
#endif

#if (TYPE_ID == DOUBLE_ID)
#pragma message "Type ==> DOUBLE"

#elif (TYPE_ID == FLOAT_ID)
#pragma message "Type ==> FLOAT"

#else
#error "You are using an undefined TYPE_ID. Only available double or float"
#endif

// Global variables
int numStepsX=0,numStepsT=0;/*!< Total number of elements in the spatial and time discretizations respectively*/
TYPE xspan=0., tspan=0., dist_mes=0;
; /*!< Total length and total time respectively*/
const TYPE pi = 3.14159265;
int *terminals; /*!< This array contains the element numbers for the boundary conditions.*/
int num_terminals; /*!< This is the number of terminals of the tree. The first element (0) is also taken into account as a terminal element.*/
TYPE **potential, *potentialR; /*!< Two arrays to save the potential.*/
TYPE *open_GL_graphE,**open_GL_graphR,***open_GL_graph; /*!< Pointers related to the graphics part.*/
char *log_file, *output; /*!< Output files*/
FILE *standard, *study; /*!< Output files*/
TYPE unit_x=0.;
int open_GL_steps=1000; /*!< Frequency to save the data*/
char *simulation;/*!< This value defines the type of simulation*/
/*! \brief This is the main function of Neurite
 */
int main(int argc, char **argv)
{
  if (argc<5 || !strcmp(argv[1],"help"))
    { 
      if(strcmp(argv[1],"help")){
	_INFO("Your number of arguments is not enough. It has to be at least 5 (type man doc/./neurite.7 )");
      }else{
	_INFO("You have some examples for running a simulation. More information you have to type man ./Neurite");
      }
      _INFO("The solver:");
      _INFO("e: explicit mode");
      _INFO("i: implicit mode");
      _INFO("Examples:");
      _INFO("%s","./Neurite e dtf 0.9 IRE 10 NRE 5 out output.txt log simulation.log");
      _INFO("%s","./Neurite i out output.txt log simulation.log");
      exit(0);
    }
  
  //Building the simulation:
  TYPE  dT=0.,time_run=0., temp=0.;
  int i=0, k=0;
  char *args;
  input in;
  params_t params = {false};
  int GL=0, GL_cont=0, interval=0;
  time_t t1, t2;
  discretization **element; // Object oriented. Class mother = discretization with two daughters = Cable and Hodgkin-Huxley
  Hodgkin_Huxley **HH_elements; // One array for each kind of elements that Neurite has
  Cable **CT_elements;
  int num_HH_elements=0;
  int num_CT_elements=0;
  TYPE E_MAX=0;
  TYPE monitoring=0;
  int dist_mes_points=0;
  
  simulation = getenv("simulation");
  if (simulation == NULL){
    _ERROR("ERROR. You did not define the specific case to simulate\n");
  }
  t1=time(NULL);

  // Initializing some variables. Default values
  numStepsX = 0; // Initial Total number of elements;
  numStepsT = 0; // Initial Total number of time steps;
  in.IRE=0;
  in.NRE =0;
  xspan =0.;
  tspan=0.;
  in.epsilon = 0.;   // Axial strain default
  in.nu=0.00261813; // Default value for an equilibrated distribution of the strain
  in.dt_factor = 1.;  //dt_factor_default
  in.init_time = 0.;
  in.E_MAX = 1;   // Maximum microscopic strain

  // handle command-line options
  args = new char[argc]();
  inputArgs(argc, argv, args, in);
  E_MAX = in.E_MAX;
  study = fopen((const char *)&output[0],"w");
  standard = fopen((const char *)&log_file[0],"w");
  args[argc-1] = '\0';
  flags_init(args, &params); // Obtaining the flags
  test_logical_options(&params,in); // testing logical options
   
  /* Here the discretization is defined. You have to define your own discretization scheme. There are some examples in the space.cpp file*/
  inputInit(in, &params, element, HH_elements, num_HH_elements, CT_elements, num_CT_elements); 
 // Strain for the electrical model
  mechanical_model_initialization(in, element, dist_mes);
  // This way of extracting the results
  dist_mes_points=measurement_point(element,dist_mes);
  // Writing some information in the log file
  fprintf(standard,"in.nu=%f\n", in.nu);
  fprintf(standard,"in.ep=%f\n", in.epsilon);
  fprintf(standard," Total time %g\n", tspan);
  fprintf(standard,"REAL The total length is =%g, dis_mes=%g\n",xspan, dist_mes);
  for(k=1;k<numStepsX-1;k++){
    fprintf(standard,"element=%i type=%i REAL epsilon =%.10f, epsilon=%g\n",k,element[k]->get_kind_HH(),in.epsilon, element[k]->get_epsilon());    
    fprintf(standard," element size = %g \n", element[k]->get_dX());
  }

  // Allocating memory
  neuriteInit(&params);
  electrical_properties(HH_elements, num_HH_elements, E_MAX);
  setting_rest_pot(potential, element);
  hhIconds(potential[0],HH_elements, num_HH_elements, E_MAX);
  cable_init_constants(CT_elements, num_CT_elements, params.HH_model);
  critical_time_step_class(HH_elements, num_HH_elements,CT_elements, num_CT_elements);
  dT = critical_time_step_selection(element);
  fprintf(standard,"stability dT=%g\n",dT);
  dT = dT*in.dt_factor;
  fprintf(standard,"applied factor %f dT=%g\n",in.dt_factor,dT);
  fprintf(standard,"total_time t span %f\n",tspan);
  temp = tspan/dT;
  numStepsT = int(temp);
  _INFO("Total time Steps=%u",numStepsT);
 
  // For open_GL
  interval = int(numStepsT/open_GL_steps);
  fprintf(standard,"intervals =%i\n", interval);
  fprintf(standard,"TotalTime=%g,Total_length=%g\n",numStepsT*dT,xspan);
  fprintf(standard,"NumStepsT=%i,NumStepsX=%i\n",numStepsT,numStepsX);
  Solver *solver;
#if defined GPU
  try {  
    _INFO("The solver is being created");
    if (params.explicit_mode){
      solver = new SolverExGPU(dT,potential,num_terminals,terminals, in.input_current,
			       num_HH_elements, HH_elements, numStepsX, element);
    }else if (params.implicit_mode) {
      solver = new SolverImGPU(dT,potential,num_terminals,terminals, in.input_current,
				 num_HH_elements, HH_elements, numStepsX, element);
    }else{
      _ERROR("You did not define the processor");
      exit(-1); 
    }

#elif defined CPU
      if (params.explicit_mode){
	solver = new SolverExCPU(dT, potential, num_terminals, terminals, 
				  in.input_current, num_HH_elements, HH_elements,
				  numStepsX, element);
      }else if (params.implicit_mode){
	solver = new SolverImCPU(dT, potential, num_terminals, terminals, 
				  in.input_current, num_HH_elements, HH_elements,
				  numStepsX, element);
      }else{
	_ERROR("ERROR: Choose explicit or implicit)");
	exit(-1); 
      }
#else
	_ERROR("You did not define the processor");
	exit(-1); 
#endif

        
      _INFO("The simulation begins...");
      int printing=0;
      int times_printing = 20; // The times you monitor the simulation in the stdout 
      monitoring=numStepsT/times_printing;

      for(i = 1; i < numStepsT-1; i++) {

	printing++;
	time_run = i * dT;

	// Updating variables
	if ((time_run >= in.init_time) && (time_run <= in.end_input)) {
	  solver->update(true);
	}else{
	  solver->update(false);
	}

	solver->calculate();

	// OUTPUT of the simulations
	// You have to be careful about how extract the results
	// Copy for open_GL
	if(++GL_cont > interval){
	  GL_cont=0;
	  solver->snapshot(open_GL_graph[0][GL++]);
	  solver->snapshot(potential[1]);
	  if(strcmp(simulation,"Ion_Channel") == 0){
	    fprintf(study,"%g %g\n",i*dT, potential[1][1]);
	  }else {
	    fprintf(study,"%.20f %.20f\n",time_run, potential[1][dist_mes_points]);
       	  }
	}
	if(printing>monitoring){
	  _INFO_REP("%.2f%% [Time = %g]",double((100*i/numStepsT)), time_run);	   
	  printing=0;
	}
      }// DT end
      _INFO_REP("%.2f%% [Time = %g]",double((100*i/numStepsT)), time_run);	   
      t2 = time(NULL);
      _INFO("\n The simulation has ended. Real_time=%f seconds", difftime(t2,t1));
      
      // Initiate and Run Glut
      solver->snapshot(potential[1]);
      if(params.openGL){
	i=0;
	unit_x = xspan/10;
	glutInit(&argc, argv);
	runGL();
      }
            
      delete[] open_GL_graphE;
      delete[] open_GL_graphR;
      delete[] open_GL_graph;
      delete[] potential;
      delete[] potentialR;
      deleting_discretization(element, HH_elements, num_HH_elements, CT_elements, num_CT_elements);
      delete[] element;
      delete[] terminals;
      delete[] HH_elements;
      delete[] CT_elements;

      delete solver;      
      delete[] args;
      fclose(standard);

#if defined GPU
    } //try
    catch (cudaException &e) {
      e.print();
    }
#endif    
    return 0;
  }

Loading data, please wait...