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;
//
//
// File author(s):  <Julian Andres Garcia Grajales>, (C) 2014
//
// Copyright: this software is licenced under the terms stipulated in the license.txt file located in its root directory
//
//


/*! 
 * \file space.cpp
 * \brief Space discretization.
 *\details This file must be edited by the user. Users give to Neurite the corresponding discretization and it is responsibility of the user to give a correct discretization filling the corresponding arrays properly. The functions shown in this file are only examples with several particular conditions.
*/
#include "space.h"
#include "configuration.h"

extern int numStepsX;
extern unsigned numStepsT;
extern double xspan, tspan, dist_mes;
extern int *terminals, num_terminals;
extern FILE *standard, *study;

int compareints2 (const void * x, const void * y)
{
  return ( *(int*)x - *(int*)y );
}

/*! Loading a neuron from the txt file. In this case, it is a passive dendritic tree.
  @param in Structure with parameters defined in the command line
  @param element Array of discretization classes
  @param CT_elements This is an array of classes that contains the cable theory elements in the discretization
  @param num_CT_elements Size of the CT_elements array

  * You have to define the neuron_for_neurite.txt file in the input folder.
 */
extern void loading_neuron_passive_dendritic_tree(input &in, discretization **&element, Cable **&CT_elements, int &num_CT_elements){

  in.input_time = 0.003; // To apply the external intensity
  in.input_current = 0.8E-9; // Amperes. External Intensity
  in.end_input = in.init_time + in.input_time;
  num_terminals = 0;
  xspan = 0; // In this case this value does not make sense
  std::vector<neuron_input> neuron;  
  // Loading the file
  int number_raws=0,i=0;
  char fname[250];
  char *execution_folder = getenv("working_directory");
  if (execution_folder != NULL){
    strcpy (fname, execution_folder);
    strcat (fname, "/neuron_for_neurite.txt");
  }else{
    fprintf(stdout,"ERROR. You did not define the working directory variable\n");
    exit(0);
  }
  //Open file
  std::ifstream inFile (fname);
 
 // Refining each element
  int refining=1; // 1-> No refinement
  // Cont terminals
  int cont_terminals=1; // This is the initial 0 element
  double micrometer=1E-6;
  
  
  if (inFile.is_open()) {
    while(!inFile.eof()){
      neuron.push_back(neuron_input());
      inFile >> neuron[i].me >> neuron[i].mother >> neuron[i].dR >> neuron[i].branching >> neuron[i].dL >> neuron[i].diameter >> neuron[i].length;
      xspan+=neuron[i].length;
      if(neuron[i].me == neuron[i].dR){
	cont_terminals++;
      }
      i++;
    } 
    number_raws = i-1;
  }else{
    // show message:
    _ERROR("Error opening file with the corresponding segmented neuron. See the space.cpp file in the loading_neuron function");
    exit(0);
  }
  num_terminals =  cont_terminals-1; 

  // Loading the neuron structure from the swc format. Pre edited with matlab to take only one branch with several branches.

 if(tspan < 1E-10){
    // Time discretization
    tspan = 0.025; // Total time
  } 

 numStepsX=number_raws*refining+1; // Total number of elements. // +1 is for the initial element
 int k=0;

  element =  new discretization*[numStepsX];
  terminals = new int[num_terminals](); // () initialize to zero
  num_CT_elements = numStepsX;
  // In this case everyone is cable theory
  num_CT_elements = numStepsX;
  CT_elements = new Cable*[num_CT_elements];
  int pos=0,cont_terminal=1;

  // Creating the fictitious 0 element
  terminals[0] = 0; // First terminal is always the fictitious element.
  pos=0;
  k++;
  CT_elements[pos]=new Cable(pos, pos + 1, -1, neuron[0].length/refining*micrometer);
  element[pos] = CT_elements[pos]; 
  element[pos]->set_epsilon(in.epsilon);
  element[pos]->set_DL(0); // No branching
  element[pos]->set_branching(false);
  element[pos]->set_diameter(neuron[0].diameter*micrometer);

 for(i=0;i<number_raws;i++){
   // If you are branching
   if(neuron[i].branching){
     for(int j=1;j<=refining;j++){
       pos = (neuron[i].me-1)*refining + j;
       CT_elements[pos]=new Cable(pos, pos + 1, neuron[i].mother*refining + j - 1, neuron[i].length/refining*micrometer);
       element[pos] = CT_elements[pos]; 
       element[pos]->set_epsilon(in.epsilon);
       if(j==refining){
	 element[pos]->set_DL(neuron[i].dL*refining - refining + 1); // branching
	 element[pos]->set_branching(true);
       }else{
	 element[pos]->set_DL(0); // branching
	 element[pos]->set_branching(false);
       }
       element[pos]->set_diameter(neuron[i].diameter*micrometer);
       k++;
     }   
   }else if(neuron[i].me == neuron[i].dR){   // If you are terminal
     for(int j=1;j<=refining;j++){
       pos = (neuron[i].me-1)*refining + j;
       if(j==refining){
       CT_elements[pos]=new Cable(pos, pos, neuron[i].mother*refining + j - 1, neuron[i].length/refining*micrometer);// The last terminal is yourself
       }else{
       CT_elements[pos]=new Cable(pos, pos+1, neuron[i].mother*refining + j - 1, neuron[i].length/refining*micrometer);
       }
       element[pos] = CT_elements[pos]; 
       element[pos]->set_epsilon(in.epsilon);
       element[pos]->set_DL(0); // No branching
       element[pos]->set_branching(false);
       element[pos]->set_diameter(neuron[i].diameter*micrometer);

       k++;
     }
     terminals[cont_terminal] = pos;
     cont_terminal++;
   }else{
     if(i>0 && neuron[i-1].me == neuron[i-1].dR){ // But the previous one was terminal
       for(int j=1;j<=refining;j++){
	 pos = (neuron[i].me-1)*refining + j;
	 if(j==1){
	   CT_elements[pos]=new Cable(pos, pos + 1, neuron[i].mother*refining + j - 1, neuron[i].length/refining*micrometer);
	 }else{
	   CT_elements[pos]=new Cable(pos, pos + 1, pos-1, neuron[i].length/refining*micrometer);
	 }
	 element[pos] = CT_elements[pos]; 
	 element[pos]->set_epsilon(in.epsilon);
	 element[pos]->set_DL(0); // No branching
	 element[pos]->set_branching(false);
	 element[pos]->set_diameter(neuron[i].diameter*micrometer);

	 k++;
       }
     }else{ // Normal normal
       for(int j=1;j<=refining;j++){
	 pos = (neuron[i].me-1)*refining + j;
	 CT_elements[pos]=new Cable(pos, pos + 1, neuron[i].mother*refining + j - 1, neuron[i].length/refining*micrometer);
	 element[pos] = CT_elements[pos]; 
	 element[pos]->set_epsilon(in.epsilon);
	 element[pos]->set_DL(0); // No branching
	 element[pos]->set_branching(false);
	 element[pos]->set_diameter(neuron[i].diameter*micrometer);
	 k++;
       }
     }
   }
 }
 dist_mes=100E-6;
 fprintf(standard,"Total steps = %i ", numStepsX);
 element[1]->set_input_current(in.input_current);
}


/*! General tree with n jumps of cable theory elements. Passive cable.
  @param in Structure with parameters defined in the command line
  @param element Array of discretization classes
  @param CT_elements This is an array of classes that contains the cable theory elements in the discretization
  @param num_CT_elements Size of the CT_elements array

  * Please, keep in your mind that you must define a new function following the way of the discretization. This function is only an example.
 */
extern void general_symmetric_tree_CT(input &in, discretization **&element, Cable **&CT_elements, int &num_CT_elements){

  TYPE dXtemp=0;
  in.input_time = 0.003; // To apply the external intensity
  in.input_current = 0.8E-9; // Nano Amperes. External Intensity
  in.init_time=0;
  in.end_input = in.init_time + in.input_time;   // when you finish to apply the current
  num_terminals = 0;
  //int min_size_branch = 70, max_size_branch = 70;// number of elements in each branch
  int *size_branch; // A pointer to know the size of the tree
  int **flag_jump, *flag_jumpE;
  int k=0;
  xspan=0;

  if(tspan < 1E-10){
    // Time discretization
    tspan = 0.025; // Total time
  } // if not is because you put as an argument.

    // All elements are CT elements. We define the element size and begin the tree in function of the elements at each branch
  dXtemp = 10E-6; //
  numStepsX=0;

  int trams=0, cont_trams=0, initial_sections = 4; // This number defines the number of jumps from the beginning.

  num_terminals = 1*pow(2,initial_sections-1) + 1;  // The last 1 is for the first fictitious element
  trams = (num_terminals-1)*2-1; // Experience

  int temp_numStepsX=0;
  size_branch = new int[trams]();// () initialize to zero
  flag_jump = new int*[trams](); // At most This number must be adjusted. THIS ALGORITHM IS TEMPORAL. flag_jump[x][y]. x is the position in the vector and y=1 if already has branch, and y=0 vice versa.
  flag_jumpE = new int[trams*2]();
   for(int i = 0; i < trams; i++){
	flag_jump[i] = flag_jumpE + (i*2);
      }

   srand(time(NULL)); // To generate random numbers. In this case, the number of elements of each branch
  // for(k=0;k<trams;k++){
  //   size_branch[k] = min_size_branch + rand()%(max_size_branch + 1 - min_size_branch);
  //   temp_numStepsX = temp_numStepsX +  size_branch[k];
  // }

   // For a fixed tree // Initial sections=4;
   size_branch[0]=44;
   size_branch[1]=62;
   size_branch[2]=48;
   size_branch[3]=50;
   size_branch[4]=28;
   size_branch[5]=50;
   size_branch[6]=26;
   size_branch[7]=32;
   size_branch[8]=48;
   size_branch[9]=50;
   size_branch[10]=30;
   size_branch[11]=32;
   size_branch[12]=40;
   size_branch[13]=55;
   size_branch[14]=38;
   temp_numStepsX = 633;

  element =  new discretization*[temp_numStepsX];
  terminals = new int[num_terminals](); // () initialize to zero

  // In this case everyone is cable theory
  num_CT_elements = temp_numStepsX ;
  CT_elements = new Cable*[num_CT_elements];

  terminals[0] = 0; // First terminal is always the fictitious element.
  k=0;
  int j=0;
  int jumps_branch = initial_sections;
  int cont_elem_tram=0, cont_jumps_tram=0, cont_back=0, cont_inside=0, cont_jumps=0;
  numStepsX=0;
  j=1;
  int mom=0;
  while(cont_trams < trams){
    cont_elem_tram = 0;
    while(cont_elem_tram < size_branch[k]){
      if(cont_elem_tram == size_branch[k]-1){
	if(cont_jumps_tram == jumps_branch-1){
	  terminals[j]=numStepsX;
	  j++;
	  CT_elements[numStepsX]=new Cable(numStepsX, numStepsX,numStepsX-1, dXtemp); //
	  element[numStepsX] = CT_elements[numStepsX]; 
	  element[numStepsX]->set_epsilon(in.epsilon);
	  cont_back = 0;
	  cont_inside=1;
	  while(cont_inside<cont_jumps){
	    cont_back++;
	    if(flag_jump[cont_jumps-1-cont_inside][1] == 0){
	      flag_jump[cont_jumps-1-cont_inside][1]=1; 
	      mom = flag_jump[cont_jumps-1-cont_inside][0];
	      element[flag_jump[cont_jumps-1-cont_inside][0]]->set_DL(numStepsX+1);
	      cont_jumps=cont_jumps-cont_inside;
	      cont_inside=10000;
	    }
	    cont_inside++;
	  }
	  cont_jumps_tram=0;
	  jumps_branch = cont_back;
	}else{
	  CT_elements[numStepsX]=new Cable(numStepsX, numStepsX+1,numStepsX-1, dXtemp);
	  element[numStepsX] = CT_elements[numStepsX]; 
	  element[numStepsX]->set_epsilon(in.epsilon);
	  element[numStepsX]->set_branching(true);
	  cont_jumps_tram++;
	}
      }else{
	if(mom == 0){
	  CT_elements[numStepsX]=new Cable(numStepsX, numStepsX+1, numStepsX-1, dXtemp);
	  element[numStepsX] = CT_elements[numStepsX]; 
	  element[numStepsX]->set_epsilon(in.epsilon);
	}else{
	  CT_elements[numStepsX]=new Cable(numStepsX, numStepsX+1, mom, dXtemp);
	  mom=0;
	  element[numStepsX] = CT_elements[numStepsX]; 
	  element[numStepsX]->set_epsilon(in.epsilon);
	}
      }
      numStepsX++;
      cont_elem_tram++;
      if(cont_elem_tram == size_branch[k]-1){
	flag_jump[cont_jumps][0]=numStepsX;
	flag_jump[cont_jumps][1]=0;
	cont_jumps++;
      }
    }
    cont_trams++;
    k++;
  }
  // Measurement distance
  dist_mes = dXtemp*numStepsX/10;

  fprintf(standard,"total elements steps=%i === %i\n", numStepsX, temp_numStepsX);
  element[1]->set_input_current(in.input_current); // Only at the initial element

}


/*! General tree with 2 jumps of cable theory elements. Passive cable. The tree is Y. This function does exactly the same than general_tree_n_terminal_CT with n=2 jumps
  @param in Structure with parameters defined in the command line
  @param element Array of discretization classes
  @param CT_elements This is an array of classes that contains the cable theory elements in the discretization
  @param num_CT_elements Size of the CT_elements array

  * Please, keep in your mind that you must define a new function following the way of the discretization. This function is only an example.
 */
extern void branching_Y_CT(input &in, discretization **&element, Cable **&CT_elements, int &num_CT_elements){

  TYPE dXtemp=0;
  in.input_time = 0.003; // To apply the external intensity
  in.input_current = 0.8E-9; // Nano Amperes. External Intensity
  in.end_input = in.init_time + in.input_time;   // when you finish to apply the current
  num_terminals =3; 
  int size_branch = 50;// number of elements for each branch.
  int temp_cont=0;
  int temp_terminals=0;

  int flag_jump=0;
  int temp_path=0; 
  xspan=0;

  if(tspan < 1E-10){
    // Time discretization
    tspan = 0.025; // Total time
  } // if not is because you put as an argument.

  
  dXtemp = 50E-6; // 50 elements for each branch with this size
  numStepsX=0;

  element =  new discretization*[size_branch*num_terminals];
  terminals = new int[num_terminals](); // () initialize to zero

  num_CT_elements = size_branch*num_terminals ;
  CT_elements = new Cable*[num_CT_elements];

  terminals[0] = 0;
  while(temp_terminals < num_terminals){
    temp_cont = 0;
    while(temp_cont<size_branch){
      temp_cont++;
      if(temp_terminals==2){
	terminals[temp_terminals-1]=numStepsX-1;

	CT_elements[numStepsX]=new Cable(numStepsX, numStepsX+1,flag_jump, dXtemp); //my mother is flag_jump
	element[numStepsX] = CT_elements[numStepsX]; 
	element[numStepsX]->set_epsilon(in.epsilon);
	element[flag_jump]->set_DL(numStepsX);
	element[flag_jump]->set_branching(true);
	element[numStepsX-1]->set_DR(numStepsX-1);
	    
	temp_terminals++;

      }else{
	CT_elements[numStepsX]=new Cable(numStepsX, numStepsX+1,numStepsX-1, dXtemp);
	element[numStepsX] = CT_elements[numStepsX]; 
	 element[numStepsX]->set_epsilon(in.epsilon);
      }
      numStepsX++;
    }
    if(temp_path==0){
      flag_jump=temp_cont-1;
    }
    temp_path++;
    temp_terminals++;
  }
  terminals[num_terminals-1] = numStepsX-1;
  element[numStepsX-1]->set_DR(numStepsX-1);
  // Measurement distance
  dist_mes = numStepsX*dXtemp/10;
  element[1]->set_input_current(in.input_current); // Only at the initial element
}

/*! General tree with 2 jumps of CT+HH elements. Active cable. The tree is Y
  @param in Structure with parameters defined in the command line
  @param element Array of discretization classes
  @param HH_elements This is an array of classes that contains the Hodgkin and Huxley elements in the discretization
  @param num_HH_elements Size of the HH_elements array
  @param CT_elements This is an array of classes that contains the cable theory elements in the discretization
  @param num_CT_elements Size of the CT_elements array

  * Please, keep in your mind that you must define a new function following the way of the discretization. This function is only an example.
 */
extern void branching_Y_CT_HH(input &in, discretization **&element, Hodgkin_Huxley **&HH_elements, int &num_HH_elements, Cable **&CT_elements, int &num_CT_elements){

  int hh_steps_temp=0, ct_steps_temp=0;
  TYPE ct_dXtemp=0, hh_dXtemp=0;
  TYPE node_size = 2.1E-6; // Node of Ranvier. From Shi papers
  TYPE inter_node = 0.0008; // internodal length. General distance
  num_terminals=3; // Initial + two final elements. Y.

  in.input_time = 0.003; // To apply the external intensity
  in.input_current = 0.8E-9; // Nano Amperes. External Intensity
  in.end_input = in.init_time + in.input_time;   // when you finish to apply the current

  if(in.IRE< 1E-10){
    // Default values for the discretization
    ct_steps_temp = 10; // internodal part
  }else{
    ct_steps_temp = in.IRE; // internodal part
  }
  if(in.NRE< 1E-10){
    // Default values for the discretization
    hh_steps_temp = 1; // internodal part
  }else{
    hh_steps_temp = in.NRE; // internodal part
  }

  if(tspan < 1E-10){
    // Time discretization
    tspan = 0.025; // Total time
  } // if not is because you put as an argument.

  //Calculate dXi, dXn and dT
  ct_dXtemp = (TYPE)inter_node/ct_steps_temp;
  hh_dXtemp = (TYPE)node_size/hh_steps_temp;

  fprintf(standard,"for discretization and simulation cabledx =%g  hhdx=%g, steps hh %i y ct %i\n", ct_dXtemp, hh_dXtemp, hh_steps_temp, ct_steps_temp);

  // Test logical options
  if(hh_dXtemp > node_size){
    _ERROR("Error. Your element size %g must be smaller than node of Ranvier size = %g\n",hh_dXtemp, node_size);
    exit(0);
  }
  if(node_size > inter_node){
    _ERROR("Error. Your node of Ranvier size = %g must be smaller than internodal distance =%g\n",node_size,inter_node);
    exit(0);
  }

  numStepsX = 0;
  fprintf(standard,"total elements steps=%i\n", numStepsX);

  int tempcts=0, temphhs=0;

  int size_branch = 50;// 
  int chunks=0;
  int temp_cont=0;
  int temp_terminals=0;


  int things=0;
  things = num_terminals*(size_branch/(hh_steps_temp+ct_steps_temp)+1);

  int flag_jump=0;
  int temp_path=0; //
  terminals = new int[num_terminals](); // () initialize to zero
  terminals[0] = 0;
  element =  new discretization*[size_branch*num_terminals+things];
  HH_elements = new Hodgkin_Huxley*[things];
  CT_elements = new Cable*[size_branch*num_terminals+things];

  while(temp_terminals < num_terminals){
    temp_cont = 0;
    while(temp_cont<size_branch){
      for(tempcts=0;tempcts<ct_steps_temp;tempcts++){
	temp_cont++;
	if(temp_terminals==2){
	  terminals[temp_terminals-1]=numStepsX-1;

	  CT_elements[num_CT_elements]=new Cable(numStepsX, numStepsX+1, flag_jump, ct_dXtemp);
	  element[numStepsX] = CT_elements[num_CT_elements];
	  num_CT_elements++;
	  element[numStepsX]->set_epsilon(in.epsilon);
	  element[flag_jump]->set_DL(numStepsX);
	  element[flag_jump]->set_branching(true);
	  element[numStepsX-1]->set_DR(numStepsX-1);
	    
	  temp_terminals++;

	}else{

	  CT_elements[num_CT_elements]=new Cable(numStepsX, numStepsX+1,numStepsX-1, ct_dXtemp);
	  element[numStepsX] = CT_elements[num_CT_elements]; 
	  num_CT_elements++;
	  element[numStepsX]->set_epsilon(in.epsilon);
	}
	numStepsX++;
      }
      for(temphhs=0;temphhs<hh_steps_temp;temphhs++){
	HH_elements[num_HH_elements] = new Hodgkin_Huxley(numStepsX, numStepsX+1,numStepsX-1, hh_dXtemp);
	element[numStepsX] = HH_elements[num_HH_elements]; 
	num_HH_elements++;
	element[numStepsX]->set_epsilon(in.epsilon);
	temp_cont++;
	numStepsX++;
      }
      chunks++;
    }
    if(temp_path==0){
      flag_jump=temp_cont-1;
    }
    temp_path++;
    temp_terminals++;
  }
   terminals[num_terminals-1] = numStepsX-1;
  element[numStepsX-1]->set_DR(numStepsX-1);
  dist_mes = numStepsX*ct_dXtemp/10;

 element[4]->set_input_current(in.input_current); // 
 element[90]->set_input_current(in.input_current); // 

}
/*! This is only for the HH model.
  @param in Structure with parameters defined in the command line
  @param element Array of discretization classes
  @param HH_elements This is an array of classes that contains the Hodgkin and Huxley elements in the discretization
  @param num_HH_elements Size of the HH_elements array
  @param CT_elements This is an array of classes that contains the cable theory elements in the discretization
  @param num_CT_elements Size of the CT_elements array

  * Please, keep in your mind that you must define a new function following the way of the discretization. This function is only an example.
 */
extern void Only_HH_Model_class(input &in, discretization **&element,  Hodgkin_Huxley **&HH_elements,int &num_HH_elements, Cable **&CT_elements, int &num_CT_elements){
  in.IRE = 1;
  in.NRE =1;
  in.dist_mes= 4.2E-3;
  xspan =6.3E-3;
  in.init_time = 0.003; // seconds. When you apply the current.
  in.input_time = 0.003;
  in.end_input = in.init_time + in.input_time;
  in.input_current = 0.8E-9;
  numStepsX =3;
  in.dt_factor = 0.0001;
  num_terminals=2;
  tspan = 0.025;
    
  element =  new discretization*[numStepsX];
  terminals = new int[num_terminals](); // () initialize to zero
  num_CT_elements = 2;
  CT_elements = new Cable*[num_CT_elements];
  num_HH_elements = 1;
  HH_elements = new Hodgkin_Huxley*[num_HH_elements];
  terminals[0] = 0;
  // First fictitious element
  CT_elements[0]=new Cable(0, 0+1, 0-1,2.1E-3);
  element[0] = CT_elements[0];
  element[0]->set_epsilon(0);
  // Active element
  HH_elements[0] = new Hodgkin_Huxley(1, 1,1-1,2.1E-3 );
  element[1] = HH_elements[0] ;
  element[1]->set_epsilon(0);
  // First fictitious element
  CT_elements[1]=new Cable(2, 1+1, 2-1,2.1E-3);
  element[2] = CT_elements[1];
  element[2]->set_epsilon(0);
  terminals[1] = 2;

 element[1]->set_input_current(in.input_current); // Only at the initial element

}

/*! General cable without branching with cable theory elements. Passive cable. This function does exactly the same than general_tree_n_terminal_CT with n=1 jumps
  @param in Structure with parameters defined in the command line
  @param element Array of discretization classes
  @param CT_elements This is an array of classes that contains the cable theory elements in the discretization
  @param num_CT_elements Size of the CT_elements array

  * Please, keep in your mind that you must define a new function following the way of the discretization. This function is only an example.
 */
extern void without_branching_CT(discretization **&element,input &in, Cable **&CT_elements, int &num_CT_elements){
  TYPE dXtemp=0;
  int cable_steps_temp=0;
  in.input_time = 0.003; // To apply the external intensity
  in.input_current = 0.8E-9; // Nano Amperes. External Intensity
  in.end_input = in.init_time + in.input_time;   // when you finish to apply the current
  num_terminals =2; // The first and last elements

  if(xspan < 1E-10){
    xspan = 0.002; // Total length. It is approximate; xspan = xspan*(1+epsilon). We put two times the total distance in order not to have boundary problems
  } // if not is because you put as an argument.

  dist_mes = xspan/10;

  if(tspan < 1E-10){
    // Time discretization
    tspan = 0.025; // Total time 
  } // if not is because you put as an argument.

  if(numStepsX != 0){
    cable_steps_temp = numStepsX;
    dXtemp = (TYPE)xspan/cable_steps_temp;
    fprintf(stdout,"dX at the whole cable =%g\n",dXtemp);
  }else{
    cable_steps_temp = 110; // Total number of elements;
    dXtemp = (TYPE)xspan/cable_steps_temp;
    numStepsX = cable_steps_temp ;
    fprintf(stdout,"dX at the whole cable =%g\n",dXtemp);
  }

  element =  new discretization*[numStepsX];
  terminals = new int[num_terminals](); // () initialize to zero
  terminals[0] = 0;
  num_CT_elements = numStepsX;
  CT_elements = new Cable*[num_CT_elements];

  for(int i=0;i<numStepsX;i++){
 
    CT_elements[i]=new Cable(i, i+1,i-1, dXtemp);
    element[i] = CT_elements[i];	  
    element[i]->set_epsilon(in.epsilon);
  }
  terminals[1] = numStepsX-1;
 element[1]->set_input_current(in.input_current); // Only at the initial element
}

/*! General cable without branching with CT+HH elements. BMMB scheme (strain, refinement, etc)
  @param in Structure with parameters defined in the command line
  @param element Array of discretization classes
  @param HH_elements This is an array of classes that contains the Hodgkin and Huxley elements in the discretization
  @param num_HH_elements Size of the HH_elements array
  @param CT_elements This is an array of classes that contains the cable theory elements in the discretization
  @param num_CT_elements Size of the CT_elements array

  * Please, keep in your mind that you must define a new function following the way of the discretization. This function is only an example.
 */
extern void without_branching_CT_HH(input &in, discretization **&element, Hodgkin_Huxley **&HH_elements,int &num_HH_elements, Cable **&CT_elements, int &num_CT_elements){

  // Temporal values
  int i=0;
  double temp=0;
  int temp2=0;
  int flag = 0, cItr=0; // to know if we are at internodal or nodal section
  int *gPoints;
  int hh_numChannels =0, cable_numCable=0, cable_numStepsX=0, hh_numStepsX=0;
  double hh_epsilon=0, cable_epsilon=0, hh_dX=0, cable_dX, cable_inter_node = 0, hh_node_size = 0;

  num_terminals =2;
  terminals = new int[num_terminals](); // () initialize to zero
  terminals[0] = 0;

  in.init_time = 0.04;
  in.input_time = 0.003; // To apply the external intensity
  in.input_current = 0.04E-9; // Nano Amperes. External Intensity.
  //in.input_current = 0;
  hh_numChannels = 0;  // to initialize the variable
  cable_numCable =0; // to initialize the variable
  in.end_input = in.init_time + in.input_time;   // when you finish to apply the current
  dist_mes = 0.010; // The referential measurement distance. BMMB paper.


  if(in.IRE< 1E-10){
    // Default values for the discretization
    //cable_numStepsX = 20; // internodal part
    cable_numStepsX = 100; // internodal part
  }else{
    cable_numStepsX = in.IRE; // internodal part
  }
  if(in.NRE< 1E-10){
    // Default values for the discretization
    hh_numStepsX = 1; // nodal part
  }else{
    hh_numStepsX = in.NRE; // nodal part
  }
  if(xspan < 1E-10){
    
    xspan = 2;//0.02 // Total length. It is approximate; xspan = xspan*(1+epsilon). We put two times the total distance in order not to have boundary problems
  } // if not is because you put as an argument.

  if(tspan < 1E-10){
    // Time discretization
    tspan = 0.040; // 0.025 Total time
  } // if not is because you put as an argument.

  if(hh_node_size < 1E-10){
    //  hh_node_size = 2.1E-6; // Node of Ranvier. From Shi papers
    hh_node_size = 2.1E-6; // Node of Ranvier. From Shi papers
  }
  if(cable_inter_node < 1E-10){
    //cable_inter_node = 0.0008; // internodal length. General distance
   cable_inter_node = 0.0008; // internodal length. General distance
  }
  //  Total number of elements
  temp2 = int((xspan/( cable_inter_node + hh_node_size))+1); // Number of CT+HH parts
  xspan = temp2*( cable_inter_node + hh_node_size); // Recalculate the total length
 
  hh_epsilon = (in.nu*xspan/(hh_node_size*temp2))*in.epsilon;
  cable_epsilon = ((1-in.nu)*xspan/(temp2*cable_inter_node))*in.epsilon;

  if(numStepsX < 1E-10){
    numStepsX = temp2*(cable_numStepsX+hh_numStepsX); // Calculate the total spatial steps
  }

  //Calculate dXi, dXn and dT
  cable_dX = (double)cable_inter_node/cable_numStepsX;
  hh_dX = (double)hh_node_size/hh_numStepsX;

  fprintf(standard,"for discretization and simulation cabledx =%g  hhdx=%g steps cable %i steps hh=%i\n", cable_dX, hh_dX, cable_numStepsX, hh_numStepsX);
  fprintf(standard,"distances cable=%.20f  hh=%.20f\n", cable_dX*cable_numStepsX, hh_dX*hh_numStepsX);

  gPoints = new int[temp2*hh_numStepsX*2](); // this 2 is to get enough memory. It does not matter, because this vector will be deleted before the simulation

  // Test logical options
  if(hh_dX > hh_node_size){
    _ERROR("Error. Your element size %g must be smaller than node of Ranvier size = %g\n",hh_dX, hh_node_size);
    exit(0);
  }
  if(hh_node_size > cable_inter_node){
    _ERROR("Error. Your node of Ranvier size = %g must be smaller than internodal distance =%g\n",hh_node_size, cable_inter_node);
    exit(0);
  }

  // be careful with this part, it is complicated
  i=1; // First ion channel is always at the first element. The element 0 is a fictitious element (BC).
  num_CT_elements++;

  temp=0;
  double tolhh=hh_dX/2;
  double tolct=cable_dX/2;
  while(i<numStepsX){
    if(flag == 0){
      gPoints[cItr] = i;
      cItr++;
      hh_numChannels++;
      temp = temp + hh_dX;
      i++;
      if(fabs(temp + tolhh) > hh_node_size){
	flag=1; // change to internodal part
	temp=0;
      }
    }
    if(flag == 1){
      i++;
      num_CT_elements++;
      temp = temp + cable_dX;
      if(fabs(temp + tolct) > cable_inter_node){
	flag=0; // change to nodal part
	temp=0;
	cable_numCable++;
      }
    } 
  }
  fprintf(standard, "Total numChannels %i. With nodal size=%f and %f as internodal distance.\n",hh_numChannels,hh_node_size, cable_inter_node);

  element =  new discretization*[numStepsX];
  num_HH_elements = hh_numChannels;
  HH_elements = new Hodgkin_Huxley*[num_HH_elements];
  CT_elements = new Cable*[num_CT_elements];
  fprintf(standard, "Total numChannels %i cable =%i total = %i numStepsX=%i\n",hh_numChannels, num_CT_elements, num_HH_elements + num_CT_elements, numStepsX);
  int tempo=0, tempoct=0;
 
  for(i=0;i<numStepsX;i++){
    if(((int*) bsearch(&i, gPoints, hh_numChannels, sizeof(int), compareints2)) != NULL){
      HH_elements[tempo] = new Hodgkin_Huxley(i, i+1,i-1, hh_dX);
      element[i] = HH_elements[tempo];
      element[i]->set_epsilon(hh_epsilon);
      tempo++; 
    }else{
      CT_elements[tempoct]=new Cable(i, i+1,i-1, cable_dX);
      element[i] = CT_elements[tempoct];
      element[i]->set_epsilon(cable_epsilon);  
      tempoct++;
    }
  }
  element[numStepsX-1]->set_DR(numStepsX-1);
  fprintf(standard,"nu=%g\n",in.nu);
  cable_numCable++;
  terminals[1] = numStepsX-1;
  fprintf(standard,"nu of equilibrium:%g\n",hh_node_size*hh_numChannels/xspan);

  fprintf(standard,"Setting the input current at the elements that you defined %g\n",in.input_current);
  element[1]->set_input_current(in.input_current); // Only at the initial element

  delete gPoints;
}




Loading data, please wait...