Increased computational accuracy in multi-compartmental cable models (Lindsay et al. 2005)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:129149
Compartmental models of dendrites are the most widely used tool for investigating their electrical behaviour. Traditional models assign a single potential to a compartment. This potential is associated with the membrane potential at the centre of the segment represented by the compartment. All input to that segment, independent of its location on the segment, is assumed to act at the centre of the segment with the potential of the compartment. By contrast, the compartmental model introduced in this article assigns a potential to each end of a segment, and takes into account the location of input to a segment on the model solution by partitioning the effect of this input between the axial currents at the proximal and distal boundaries of segments. For a given neuron, the new and traditional approaches to compartmental modelling use the same number of locations at which the membrane potential is to be determined, and lead to ordinary differential equations that are structurally identical. However, the solution achieved by the new approach gives an order of magnitude better accuracy and precision than that achieved by the latter in the presence of point process input.
Reference:
1 . Lindsay AE, Lindsay KA, Rosenberg JR (2005) Increased computational accuracy in multi-compartmental cable models by a novel approach for precise point process localization. J Comput Neurosci 19:21-38 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s):
Channel(s): I Na,t; I K;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; C or C++ program;
Model Concept(s): Methods;
Implementer(s):
Search NeuronDB for information about:  I Na,t; I K;
/
LindsayEtAl2005
readme.txt
03-192.pdf
AnalyseResults.c
BitsAndPieces.c
CellData.dat
CompareSpikeTrain.c
Ed04.tex
ExactSolution.dat
GammaCode
Gen.tex
Gen1.tex
Gen2.tex
Gen3.tex
Gen4.tex
Gen5.tex
Gen6.tex
GenCom.c
GenCom1.c
GenCom2.c
GenComExactSoln.c
GenerateInput.c
GenerateInputText.c
GenRan.ran
GetNodeNumbers.c
Info100.dat
Info20.dat
Info200.dat
Info30.dat
Info300.dat
Info40.dat
Info400.dat
Info50.dat
Info500.dat
Info60.dat
Info70.dat
Info80.dat
Info90.dat
InputCurrents.dat
InputDendrite.dat
JaySpikeTrain.c
JayTest1.dat
JayTest100.dat
KenSpikeTrain.c
KenTest1.dat *
KenTest10.dat
KenTest100.dat *
KenTest10p.dat
KenTest1p.dat *
KenTest2.dat
KenTest2p.dat
KenTest3.dat
KenTest3p.dat
KenTest4.dat
KenTest4p.dat
KenTest5.dat
KenTest5p.dat
KenTest6.dat
KenTest6p.dat
KenTest7.dat
KenTest7p.dat
KenTest8.dat
KenTest8p.dat
KenTest9.dat
KenTest9p.dat
LU.c
Mean50.dat
Mean500.dat
mosinit.hoc
NC.pdf
NC.tex
NC1.tex
NC2.tex
NC3.tex
NC4.tex
NC5.tex
NC6.tex
NCFig2.eps *
NCFig3.eps *
NCFig4.eps *
NCFig5a.eps *
NCFig5b.eps *
NCFig6.eps *
NCPics.tex
NeuronDriver.hoc
NewComExactSoln.c
NewComp.pdf
NewComp.ps
NewComp.tex
NewComp.toc
NewComp1.tex
NewComp2.tex
NewComp3.tex
NewComp4.tex
NewComp5.tex
NewComp6.tex
NewCompFig1.eps
NewCompFig2.eps *
NewCompFig3.eps *
NewCompFig4.eps *
NewCompFig5a.eps *
NewCompFig5b.eps *
NewCompFig6.eps *
NewCompPics.tex
NewComSpikeTrain.c
NewRes.dat
NewRes60.dat
NewRes70.dat
NewRes80.dat
NewSynRes40.dat
NewTestCell.d3
NResults.res
OldComExactSoln.c
out.res
principles_01.tex
rand
Ratio.dat
RelErr.dat
ReviewOfSpines.pdf
SpikeTimes.dat
TestCell.d3
TestCell1.d3
TestCell2.d3
TestCell3.d3
TestCell4.d3
testcellnew2.hoc
TestCGS.c
TestGen1.c
TestSim.hoc
TestSim020.hoc
TestSim030.hoc
TestSim040.hoc
TestSim050.hoc
TestSim060.hoc
TestSim070.hoc
TestSim080.hoc
TestSim090.hoc
TestSim1.hoc
TestSim100.hoc
TestSim200.hoc
TestSim300.hoc
TestSim400.hoc
TestSim500
TestSim500.hoc
                            
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

 /***************************************************************
    Function to analyse the numerical error of Generalised
    Compartmental Model. A single input is simulated and
    the exact solution determined using the Equivalent Cable
  ***************************************************************/

/* Global definitions */
#define             CS    1.0
#define             GS    0.091
#define             GA    14.286
#define             CM    1.0
#define             GM    0.091
#define         OUTPUT    "KenSpike200.dat"
#define           TEND    11        /* Must be an integer */
#define            TAU    0.5
#define            CGS    1.0e-26   /* Tolerance used in CGS algorithm */
#define           GMAX    3.0e-5
#define           RATE    30.0
#define           VSYN    115.0
#define             NT    1000
#define          NODES    200
#define          NSEED    2         /* Seed for random number generator */
#define        CELSIUS    18.5      /* Celsius temperature of neuron */

/* Parameters for exact solution */
#define           NCON    100        /* Number of contacts */
#define             RS    0.005

typedef struct SparseMatrix_t
{
        double *a;
        int *col;
        int *StartRow;
        int n;

} SparseMatrix;


typedef struct cond_t {
        int n;               /* No. timesteps in conductance profile */
        double dt;           /* Integration time step (msecs) */
        double *g;           /* Conductance profile for prescribed dt (mu a) */
        double tol;          /* Determines off-threshold for synapse */
        double tau;          /* Rise time for synapse (msecs) */
        double gmax;         /* Peak conductance per synapse (mu) */
} cond;


typedef struct synapse_t
{
        int id;                     /* Identifies contact type */
        double vlam;                /* Location of contact */

        int np;                     /* Proximal neighbour */
        int nd;                     /* Distal neighbour */

/*  Properties of synaptic conductance profile */
        cond *SynCond;              /* Address of synaptic conductance profile */
        int SynTime;                /* Time to activation of synapse/time since activation */
        int MaxStep;                /* Time of inactivation of synapse */
        double vsyn;                /* Reversal potential of synaptic species */
        struct synapse_t *NextSyn;  /* Address of next synapse */
} synapse;


typedef struct segment_t
{
/*  Physical properties of segment */
        double rp;                  /* Proximal radius of segment */
        double rd;                  /* Distal radius of segment */
        double hseg;                /* Dendritic segment length (cm) */
        int nsyn;                   /* Number of synapses on segment */

/*  Synaptic conductances */
        double *gold;               /* Old conductance */
        double *gnew;               /* New conductance */

/*  Auxiliary vectors */
        double *vlam;
        double *diag;
        double *vec;
        double *phi;
        double *phi1;
        double *phi2;
        double *phi3;

/*  Solution vector */
        double *soln;

/*  Auxiliary scalars */
        double OldProx1;
        double OldProx2;
        double OldProx3;
        double NewProx1;
        double NewProx2;
        double NewProx3;
        double OldDist1;
        double OldDist2;
        double OldDist3;
        double NewDist1;
        double NewDist2;
        double NewDist3;

/*  Contact information */
        synapse *SynList;           /* List of segment synapses */
} segment;



typedef struct branch_t
{
/*  Connectivity of branch */
        struct branch_t *parent;    /* Address of parent branch */
        struct branch_t *child;     /* Address of child branch */
        struct branch_t *peer;      /* Addresss of peer branch */

/*  Physical properties of branch */
        int id;                     /* Branch identifier */
        int nseg;                   /* Number of segments specifying branch */
        double xl;                  /* X-coordinate of lefthand endpoint */
        double yl;                  /* Y-coordinate of lefthand endpoint */
        double zl;                  /* Z-coordinate of lefthand endpoint */
        double xr;                  /* X-coordinate of righthand endpoint */
        double yr;                  /* Y-coordinate of righthand endpoint */
        double zr;                  /* Z-coordinate of righthand endpoint */
        double diam;                /* Branch diameter (cm) */
        double plen;                /* Branch length (cm) */
        double hseg;                /* Dendritic segment length (cm) */

/*  Node information for spatial representation */
        int nodes;                  /* Total number nodes spanning branch */
        int junct;                  /* Junction node of the branch */
        int first;                  /* Internal node connected to junction */

/*  Contact information */
        segment *SegList;           /* Pointer to an array of segments */
} branch;


typedef struct dendrite_t
{
        branch *root;               /* Pointer to root branch of dendrite */
        double plen;                /* Length of dendrite */
} dendrite;


typedef struct neuron_t
{
        int ndend;                  /* Number of dendrites */
        dendrite *dendlist;         /* Pointer to an array of dendrites */
} neuron;


/* Function type declarations */
cond    *ConductanceProfile( double,  double, double, double );

int     Count_Synapses( branch *, branch *);

double  branch_length( branch *, branch *),
        ran(unsigned int *, unsigned int *, unsigned int *),
        alfa_h( double ),
        alfa_m( double ),
        alfa_n( double ),
        beta_h( double ),
        beta_m( double ),
        beta_n( double );

void    Build_Test_Dendrite( branch **, branch *),
        Output_Information( branch *, FILE *),
        Remove_Branch( branch **, branch *),
        Enumerate_Nodes( branch *, int *),
        Generate_Dendrite(branch *, int *),
        Initialise_Synapses( branch *),
        Update_Synapses( branch *),
        solve( int, double *, double *, double *),
        Matrix_Vector_Multiply( SparseMatrix *, double *, double *),
        Matrix_Malloc( SparseMatrix *, int, int),
        Matrix_Free( SparseMatrix *),
        Update_Sparse_Matrices( branch *, int *),
        cgs( int *, SparseMatrix *, double *, double *, double);


/* Global Variables */
SparseMatrix lhs, rhs;
double pi, dt, *NewAmp;
unsigned int ix, iy, iz;

int main( int argc, char **argv )
{
    extern unsigned int ix, iy, iz;
    extern SparseMatrix lhs, rhs;
    extern double pi, dt, *NewAmp;
    int k, j, id, nn, nodes, n, nseg, i, in, nstep, maxstep, FirstNode,
        nsyn, NumberOfBranches, counter, connected, spk, nspk, getmem;
    double *v, *x, *OldAmp, max, AreaOfSoma, sum, tmp, vs, len, h, sc,
           gs, interval, dx, CellLength, LocusContact, mval, nval,
           hval, aval, bval, vm0, vm1, vm2, tnow;
    double *StoredLHS, *StoredRHS;
    double v_na=115.0, v_k=-12.0, v_l, g_na=120.0, g_k=36.0, g_l=0.3;
    void srand( unsigned int);
    neuron *cell;
    cond *SynCond;
    synapse *NewSyn, *syn;
    segment *seg;
    branch *bo, *bn, *FirstBranch;
    char word[20];
    FILE *fp, *fp1;

/*  Initialise random number generator */
    fp = fopen(OUTPUT,"w");
    fclose(fp);
    nspk = spk = 0;
    dt = 1.0/((double) NT);
    srand( ((unsigned int) NSEED) );
    ix = rand( );
    iy = rand( );
    iz = rand( );
    SynCond = ConductanceProfile( dt, TAU, 5.e-7, GMAX );

/*  Load Test Neuron */
    maxstep =  1000*NT*TEND;
    pi = 4.0*atan(1.0);
    AreaOfSoma = 4.0*pi*RS*RS;
    if ( argc != 2 ) {
        printf("\n Invoke program with load <input>\n");
        return 1;
    } else {
        printf("\nOpening file %s\n",argv[1]);
        if ( (fp=fopen(argv[1],"r")) == NULL ) {
            printf("\n Test Neuron file not found");
            return 1;
        }
    }

/*  Get branch data */
    bo = NULL;
    while  ( fscanf(fp,"%s",word) != EOF ) {
        if ( strcmp(word,"Branch") == 0 || strcmp(word,"branch") == 0 ) {
            bn = (branch *) malloc( sizeof(branch) );
            fscanf(fp,"%d", &(bn->id) );
            bn->peer = NULL;
            bn->child = NULL;
            if ( bo != NULL) {
                bo->child = bn;
            } else {
                FirstBranch = bn;
            }
            bn->parent = bo;
            fscanf(fp,"%lf %lf %lf", &(bn->xl), &(bn->yl), &(bn->zl) );
            fscanf(fp,"%lf %lf %lf", &(bn->xr), &(bn->yr), &(bn->zr) );
            fscanf(fp,"%lf %lf", &(bn->plen), &(bn->diam) );
            bo = bn;
        } else {
            printf("Unrecognised dendritic feature\n");
            return 0;
        }
    }
    fclose(fp);

/*  STEP 1A. - Determine cell length and number of branches */
    CellLength = 0.0;
    NumberOfBranches = 0;
    bn = FirstBranch;
    while ( bn ) {
        NumberOfBranches++;
        CellLength += bn->plen;
        bn = bn->child;
    }

/*  STEP 1B. - Determine segment lengths and allocate segments */
    h = CellLength/((double) NODES-NumberOfBranches);
    bn = FirstBranch;
    while ( bn ) {
        nseg = bn->nseg = ((int) ceil((bn->plen)/h));
        bn->hseg = (bn->plen)/((double) bn->nseg);
        bn->SegList = (segment *) malloc( nseg*sizeof(segment) );
        for ( k=0 ; k<nseg ; k++ ) {
            (bn->SegList)[k].hseg = bn->hseg;
            (bn->SegList)[k].rp = 0.5*(bn->diam);
            (bn->SegList)[k].rd = 0.5*(bn->diam);
            (bn->SegList)[k].SynList = NULL;
            (bn->SegList)[k].gold = NULL;
            (bn->SegList)[k].gnew = NULL;
            (bn->SegList)[k].soln = NULL;
            (bn->SegList)[k].nsyn = 0;
        }
        bn = bn->child;
    }

/*  STEP 1C. - Randomly place NCON synapses on branches */
    for ( nsyn=0 ; nsyn<NCON ; nsyn++ ) {
        LocusContact = CellLength*ran( &ix, &iy, &iz);
        bn = FirstBranch;
        len = bn->plen;
        while ( LocusContact > len ) {
            bn = bn->child;
            len += bn->plen;
        }
        NewSyn = (synapse *) malloc( sizeof(synapse) );
        NewSyn->NextSyn = NULL;
        NewSyn->SynCond = SynCond;
        NewSyn->MaxStep = SynCond->n;
        NewSyn->vsyn = VSYN;
        LocusContact -= (len-(bn->plen));
        LocusContact /= bn->hseg;
        NewSyn->vlam = fmod(LocusContact,1.0);
        nseg = ((int) floor(LocusContact));
        (bn->SegList)[nseg].nsyn++;
        syn = (bn->SegList)[nseg].SynList;
        if ( syn ) {
            if ( NewSyn->vlam < syn->vlam ) {
                NewSyn->NextSyn = syn;
                (bn->SegList)[nseg].SynList = NewSyn;
            } else {
                while ( syn->NextSyn && NewSyn->vlam > syn->NextSyn->vlam ) syn = syn->NextSyn;
                if ( syn->NextSyn ) NewSyn->NextSyn = syn->NextSyn;
                syn->NextSyn = NewSyn;
            }
        } else {
            (bn->SegList)[nseg].SynList = NewSyn;
        }
    }

/*  STEP 2A. - Count root branches */
    bo = FirstBranch;
    n = 0;
    while ( bo ) {
        bn = FirstBranch;
        do {
            tmp = pow(bo->xl-bn->xr,2)+pow(bo->yl-bn->yr,2)+
                  pow(bo->zl-bn->zr,2);
            connected = ( tmp < 0.01 );
            bn = bn->child;
        } while ( bn && !connected );
        if ( !connected ) n++;
        bo = bo->child;
    }

/*  STEP 2B. - Identify somal dendrites but extract nothing */
    printf("\nTree contains %d individual dendrite(s) ...\n", n);
    cell = (neuron *) malloc( sizeof(neuron) );
    cell->ndend = n;
    cell->dendlist = (dendrite *) malloc( n*sizeof(dendrite) );
    bo = FirstBranch;
    n = 0;
    while ( n < cell->ndend ) {
        bn = FirstBranch;
        do {
            tmp = pow(bo->xl-bn->xr,2)+pow(bo->yl-bn->yr,2)+
                  pow(bo->zl-bn->zr,2);
            connected = ( tmp < 0.01 );
            bn = bn->child;
        } while ( bn );
        if ( !connected ) cell->dendlist[n++].root = bo;
        bo = bo->child;
    }

/*  STEP 2C. - Extract root of each dendrite from dendrite list */
    for ( k=0 ; k<cell->ndend ; k++ ) {
        bo = cell->dendlist[k].root;
        Remove_Branch( &FirstBranch, bo);
    }

/*  STEP 2D. - Build each test dendrite from its root branch */
    for ( k=0 ; k<cell->ndend ; k++ ) {
        Build_Test_Dendrite( &FirstBranch, cell->dendlist[k].root );
    }
    if ( FirstBranch != NULL ) printf("\nWarning: Unconnected branch segments still exist\n");

/*  STEP 2E. - Count number of synapses on Cell */
    for ( nsyn=k=0 ; k<cell->ndend ; k++ ) {
        bn = cell->dendlist[k].root;
        nsyn += Count_Synapses( cell->dendlist[k].root, bn);
    }
    printf("\nNumber of Synapses %d", nsyn);

/*  STEP 3A. - Enumerate Nodes */
    FirstNode = 0;
    for ( k=0 ; k<cell->ndend ; k++ ) Enumerate_Nodes( cell->dendlist[k].root, &FirstNode );
    for ( k=0 ; k<cell->ndend ; k++ ) {
        cell->dendlist[k].root->junct = FirstNode;
    }
    printf("\nNumber of nodes is %d\n", FirstNode+1);
    getchar( );

/*  STEP 3B. - Construct Sparse Matrices */
    nodes = FirstNode+1;
    Matrix_Malloc( &lhs, nodes, 3*nodes-2 );
    Matrix_Malloc( &rhs, nodes, 3*nodes-2 );
    lhs.StartRow[0] = rhs.StartRow[0] = 0;
    for ( counter=k=0 ; k<cell->ndend ; k++ ) {
        bn = cell->dendlist[k].root;
        Generate_Dendrite( bn, &counter);
    }
    lhs.n = rhs.n = nodes;

/*  STEP 3C. - Do somal node */
    lhs.a[3*nodes-3] = rhs.a[3*nodes-3] = 0.0;
    for ( k=0 ; k<cell->ndend ; k++ ) {
        bn = cell->dendlist[k].root;
        lhs.a[counter] = (bn->diam)*(bn->hseg)/6.0;
        rhs.a[counter] = -0.25*pi*pow(bn->diam,2)/(bn->hseg);
        lhs.col[counter] = rhs.col[counter] = bn->first;
        lhs.a[3*nodes-3] += 2.0*pi*(bn->diam)*(bn->hseg)/6.0;
        rhs.a[3*nodes-3] += 0.25*pi*pow(bn->diam,2)/(bn->hseg);
        counter++;
    }
    lhs.col[counter] = rhs.col[counter] = nodes-1;
    lhs.StartRow[nodes] = rhs.StartRow[nodes] = counter+1;

/*  STEP 3D. - Construct Vectors to hold currents */
    OldAmp = (double *) malloc( nodes*sizeof(double) );
    NewAmp = (double *) malloc( nodes*sizeof(double) );
    for ( k=0 ; k<nodes ; k++ ) NewAmp[k] = 0.0;

/*  STEP 3E. - Fill in properties of synapses */
    for( k=0 ; k<cell->ndend ; k++ ) {
        bn = cell->dendlist[k].root;
        Initialise_Synapses(bn);
    }
    for ( k=0 ; k<3*nodes-2 ; k++ ) {
        rhs.a[k] = 0.5*dt*(GA*rhs.a[k]+GM*lhs.a[k]);
        rhs.a[k] = CM*lhs.a[k]-rhs.a[k];
        lhs.a[k] = 2.0*CM*lhs.a[k]-rhs.a[k];
    }
    lhs.a[3*nodes-3] += AreaOfSoma*CS;
    rhs.a[3*nodes-3] += AreaOfSoma*CS;

/*  STEP 13. - Construct and load vectors to hold transient information */
    StoredLHS = (double *) malloc( (3*nodes-2)*sizeof(double) );
    StoredRHS = (double *) malloc( (3*nodes-2)*sizeof(double) );
    for ( k=0 ; k<3*nodes-2 ; k++ ) {
        StoredLHS[k] = lhs.a[k];
        StoredRHS[k] = rhs.a[k];
    }

//    fp1 = fopen("KenOutput.dat","w");
//    for( k=0 ; k<3*nodes-2 ; k++ ) {
//        fprintf(fp1,"%12.10lf\n", lhs.a[k]);
//    }
//    fclose(fp1);
//    printf("\nOutput information");
//    getchar( );

/*  STEP 14. - Compute somal conductances and the leakage potential */
    g_na *= AreaOfSoma;
    g_k *= AreaOfSoma;
    g_l *= AreaOfSoma;
    hval = alfa_h(0.0)/(alfa_h(0.0)+beta_h(0.0));
    mval = alfa_m(0.0)/(alfa_m(0.0)+beta_m(0.0));
    nval = alfa_n(0.0)/(alfa_n(0.0)+beta_n(0.0));
    v_l = g_na*pow(mval,3)*hval*v_na+g_k*pow(nval,4)*v_k;
    v_l = -v_l/g_l;

/*  STEP 15. - Construct and initialise potentials and currents */
    v = (double *) malloc( (nodes)*sizeof(double) );
    x = (double *) malloc( (nodes)*sizeof(double) );
    for ( k=0 ; k<nodes ; k++ ) v[k] = 0.0;

/*  Initialise temporal integration and integrate forward */
    nstep = 0;
    getmem = 1;
    while ( nstep < maxstep ) {
        nstep++;
//        if ( nstep%1000 == 0 ) printf("\r Msec %d",nstep/1000);

/*  Phase 1. - Update HH channel variables */
        vs = v[nodes-1];
        aval = dt*alfa_h(vs);
        bval = dt*beta_h(vs);
        tmp = 0.5*(aval+bval);
        hval = (aval+(1.0-tmp)*hval)/(1.0+tmp);
        aval = dt*alfa_m(vs);
        bval = dt*beta_m(vs);
        tmp = 0.5*(aval+bval);
        mval = (aval+(1.0-tmp)*mval)/(1.0+tmp);
        aval = dt*alfa_n(vs);
        bval = dt*beta_n(vs);
        tmp = 0.5*(aval+bval);
        nval = (aval+(1.0-tmp)*nval)/(1.0+tmp);

/*  Phase 2. - Compute somal conductance and contribution to current */
        gs = g_l;                       /* Leakage conductance */
        sc = g_l*v_l;                   /* Leakage contribution to somal current */
        tmp = g_na*hval*pow(mval,3);    /* Sodium conductance */
        gs += tmp;
        sc += tmp*v_na;                 /* Sodium contribution to somal current */
        tmp = g_k*pow(nval,4);          /* Potassium conductance */
        gs += tmp;
        sc += tmp*v_k;                  /* Potasium contribution to somal current */

/*  Phase 3. - Zero LHS, RHS and synaptic currents */
        for ( k=0 ; k<3*nodes-2 ; k++ ) lhs.a[k] = rhs.a[k] = 0.0;
        for ( k=0 ; k<nodes ; k++ ) {
            OldAmp[k] = NewAmp[k];
            NewAmp[k] = 0.0;
        }

/*  Phase 4. - Update synaptic conductances and input */
        counter = 0;
        for ( k=0 ; k<cell->ndend ; k++ ) {
            bn = cell->dendlist[k].root;
            Update_Synapses( bn );
            bn = cell->dendlist[k].root;
            Update_Sparse_Matrices( bn, &counter);
        }

/*  Phase 4a. - Update soma */
        tmp = 0.5*dt;
        for ( k=0 ; k<cell->ndend ; k++ ) {
            bn = cell->dendlist[k].root;
            seg = &(bn->SegList)[0];
            lhs.a[counter] = tmp*(seg->NewProx2);
            rhs.a[counter] = -tmp*(seg->OldProx2);
            NewAmp[nodes-1] += seg->NewProx3;
            lhs.a[3*nodes-3] += tmp*(seg->NewProx1);
            rhs.a[3*nodes-3] -= tmp*(seg->OldProx1);
            counter++;
        }

/*  Phase 5. - Complete the construction of LHS and RHS matrices */
        for ( k=0 ; k<3*nodes-2 ; k++ ) {
            lhs.a[k] += StoredLHS[k];
            rhs.a[k] += StoredRHS[k];
        }
        gs *= 0.5*dt;
        lhs.a[3*nodes-3] += gs;
        rhs.a[3*nodes-3] -= gs;

/*  Phase 6. - Step potential forward */
        Matrix_Vector_Multiply( &rhs, v, x);
        x[nodes-1] += sc*dt;
        tmp = 0.5*dt;
        for ( k=0 ; k<nodes ; k++ ) x[k] -= tmp*(NewAmp[k]+OldAmp[k]);
        cgs( &getmem, &lhs, x, v, CGS);

/*  Phase 7. - Test for spikes */
        if ( nstep == 1 ) {
            vm2 = 0.0;
            vm1 = v[nodes-1];
        } else {
            vm0 = v[nodes-1];
            if ( !spk ) {
                spk = ( vm0 > 50.0 && vm1 > vm2 && vm1 > vm0 );
                if ( spk ) {
                    if ( nstep >= 1000*NT ) {
                        nspk++;
                        tnow = dt*((double) nstep)-1000.0;
                        tmp = tnow+0.5*dt*(vm2+3.0*vm0-4.0*vm1)/(vm0-2.0*vm1+vm2);
                        nn = ((int) floor(tmp));
                        if ( fmod(tmp,1.0)>0.5 ) nn++;
                        fp = fopen(OUTPUT,"a");
                        fprintf(fp,"%d\n",nn);
                        fclose(fp);
                    } else {
                        printf("\nUnrecorded spike fired\n");
                    }
                }
            }
            vm2 = vm1;
            vm1 = vm0;
        }

/*  Phase 8. - Reset spike flag */
        if ( vs < 0.0 && spk == 1 ) spk = 0;
        if ( nstep%(500*NT) == 0 ) {
            tnow = dt*((double) nstep/NT);
            printf("\rReached time %5.1lf ms \t Spikes so far %d", tnow, nspk);
        }
    }
    return 0;
}


 /******************************************************
     Function to build a test dendrite from its root
  ******************************************************/
void Build_Test_Dendrite( branch **head, branch *root)
{
    double tmp;
    branch *bnow, *bnext, *btmp;

    bnow = *head;
    while ( bnow != NULL ) {

/*  Store bnow's child in case it's corrupted */
        bnext = bnow->child;

/*  Decide if proximal end of bnow is connected to distal end of root */
        tmp = pow(bnow->xl-root->xr,2)+
              pow(bnow->yl-root->yr,2)+
              pow(bnow->zl-root->zr,2);
        if ( tmp <= 0.01 ) {

/*  Remove bnow from the branch list */
            Remove_Branch( head, bnow);

/*  Connect bnow to the root as the child or a peer of the child.
    Initialise childs' children and peers to NULL as default */
            bnow->child = NULL;
            bnow->peer = NULL;
            bnow->parent = root;

/*  Inform root about its child if it's the first child, or add
    new child to first child's peer list */
            if ( root->child != NULL ) {
                btmp = root->child;
                while ( btmp->peer != NULL ) btmp = btmp->peer;
                btmp->peer = bnow;
            } else {
                root->child = bnow;
            }
        }

/*  Initialise bnow to next branch in list */
        bnow = bnext;
    }

/* Iterate through remaining tree */
    if ( root->child ) Build_Test_Dendrite( head, root->child);
    if ( root->peer ) Build_Test_Dendrite( head, root->peer);
    return;
}


 /*********************************************************
        Function to remove a branch from a branch list
  *********************************************************/
void Remove_Branch(branch **head, branch *b)
{
    if ( *head == NULL || b == NULL ) return;
    if ( *head == b ) {
        *head = b->child;
        if ( *head != NULL )  (*head)->parent = NULL;
    } else {
        b->parent->child = b->child;
        if ( b->child != NULL ) b->child->parent = b->parent;
    }
    b->parent = NULL;
    b->child = NULL;
    return;
}



 /*********************************************
      Function to count synapses on a branch
  *********************************************/
int Count_Synapses( branch *bstart, branch *bnow)
{
    static int n;
    int k;
    synapse *syn;

    if ( bstart == bnow ) n = 0;
    if ( bnow->child ) Count_Synapses(bstart, bnow->child);
    if ( bnow->peer ) Count_Synapses(bstart, bnow->peer);

    for ( k=0 ; k<bnow->nseg ; k++ ) {
        syn = (bnow->SegList)[k].SynList;
        while ( syn ) {
            n++;
            syn = syn->NextSyn;
        }
    }
    return n;
}


 /*******************************************************
       Function to enumerate the nodes on a dendrite
  *******************************************************/
void Enumerate_Nodes(branch *bnow, int *FirstNode )
{
    int node, k;
    branch *btmp;

    if ( bnow->child ) Enumerate_Nodes( bnow->child, FirstNode );
    if ( bnow->peer ) Enumerate_Nodes( bnow->peer, FirstNode );

    if ( bnow->child ) {
        btmp = bnow->child;
        while ( btmp ) {
            btmp->junct = *FirstNode;
            btmp = btmp->peer;
        }
    }
    *FirstNode += bnow->nseg;
    bnow->first = *FirstNode-1;
    return;
}


 /*******************************************************
                 Function to output information
  *******************************************************/
void Output_Information(branch *b, FILE *fp)
{
    int k;
    double pos;
    synapse *syn;
    segment *seg;

    if ( b->child ) Output_Information( b->child, fp);
    if ( b->peer ) Output_Information( b->peer, fp);

    for ( k=0 ; k<b->nseg ; k++ ) {
        seg = &(b->SegList)[k];
        syn = seg->SynList;
        while ( syn ) {
            pos = syn->vlam+((double) k);
            fprintf(fp," %12.6lf \t %d \n",pos, syn->SynTime);
            syn = syn->NextSyn;
        }
    }
    fprintf(fp,"End of branch\n");
    return;
}


 /***************************************************
        Function to constuct sparse matrices
  ***************************************************/
void Generate_Dendrite( branch *b, int *counter)
{
    int k, CurrentNode, nc;
    extern double pi;
    extern SparseMatrix lhs, rhs;
    branch *btmp;
    double SumL, SumR;

/* Step 1 - Recurse to the end of the dendrite */
    if ( b->child ) Generate_Dendrite( b->child, counter);
    if ( b->peer ) Generate_Dendrite( b->peer, counter);

/* Step 2 - Build matrix entries for distal node of branch */
    nc = b->nseg;
    CurrentNode = (b->first)-(nc-1);
    if ( b->child ) {
        btmp = b->child;
        SumR = SumL = 0.0;
        while ( btmp ) {
            lhs.a[*counter] = pi*(btmp->diam)*(btmp->hseg)/6.0;
            rhs.a[*counter] = -0.25*pi*pow(btmp->diam,2)/(btmp->hseg);
            SumL += 2.0*pi*(btmp->diam)*(btmp->hseg)/6.0;
            SumR += 0.25*pi*pow(btmp->diam,2)/(btmp->hseg);
            lhs.col[*counter] = rhs.col[*counter] = btmp->first;
            (*counter)++;
            btmp = btmp->peer;
        }
        lhs.a[*counter] = SumL+2.0*pi*(b->diam)*(b->hseg)/6.0;
        rhs.a[*counter] = SumR+0.25*pi*pow(b->diam,2)/(b->hseg);
        lhs.col[*counter] = rhs.col[*counter] = CurrentNode;
        (*counter)++;
        lhs.a[*counter] = pi*(b->diam)*(b->hseg)/6.0;
        rhs.a[*counter] = -0.25*pi*pow(b->diam,2)/(b->hseg);
        if ( CurrentNode == b->first ) {
            lhs.col[*counter] = rhs.col[*counter] = b->junct;
        } else {
            lhs.col[*counter] = rhs.col[*counter] = CurrentNode+1;
        }
        (*counter)++;
        lhs.StartRow[CurrentNode+1] = rhs.StartRow[CurrentNode+1] = *counter;
    } else {
        lhs.a[*counter] = 2.0*pi*(b->diam)*(b->hseg)/6.0;
        rhs.a[*counter] = 0.25*pi*pow(b->diam,2)/(b->hseg);
        lhs.col[*counter] = rhs.col[*counter] = CurrentNode;
        (*counter)++;
        lhs.a[*counter] = pi*(b->diam)*(b->hseg)/6.0;
        rhs.a[*counter] = -0.25*pi*pow(b->diam,2)/(b->hseg);
        if ( CurrentNode == b->first ) {
            lhs.col[*counter] = rhs.col[*counter] = b->junct;
        } else {
            lhs.col[*counter] = rhs.col[*counter] = CurrentNode+1;
        }
        (*counter)++;
        lhs.StartRow[CurrentNode+1] = rhs.StartRow[CurrentNode+1] = *counter;
    }

/* Step 3 - Build matrix entries for internal nodes of branch */
    for ( k=nc-1 ; k>0 ; k-- ) {
        CurrentNode++;
        lhs.a[*counter] = pi*(b->diam)*(b->hseg)/6.0;
        rhs.a[*counter] = -0.25*pi*pow(b->diam,2)/(b->hseg);
        lhs.col[*counter] = rhs.col[*counter] = CurrentNode-1;
        (*counter)++;
        lhs.a[*counter] = 4.0*pi*(b->diam)*(b->hseg)/6.0;
        rhs.a[*counter] = 0.5*pi*pow(b->diam,2)/(b->hseg);
        lhs.col[*counter] = rhs.col[*counter] = CurrentNode;
        (*counter)++;
        lhs.a[*counter] = pi*(b->diam)*(b->hseg)/6.0;
        rhs.a[*counter] = -0.25*pi*pow(b->diam,2)/(b->hseg);
        if ( CurrentNode == b->first ) {
            lhs.col[*counter] = rhs.col[*counter] = b->junct;
        } else {
            lhs.col[*counter] = rhs.col[*counter] = CurrentNode+1;
        }
        (*counter)++;
        lhs.StartRow[CurrentNode+1] = rhs.StartRow[CurrentNode+1] = *counter;
    }
    return;
}


 /***********************************************
       Function to assign synaptic weights
  ***********************************************/
void Initialise_Synapses( branch *b )
{
    extern double pi;
    int j, jj, k, nsyn;
    double fac, interval, vold, vnew;
    synapse *syn;
    segment *seg;

    if ( b->child ) Initialise_Synapses( b->child );
    if ( b->peer ) Initialise_Synapses( b->peer );

    for ( k=0 ; k<b->nseg ; k++ ) {

/*  Phase 1. - Allocate storage for synaptic activity */
        nsyn = (b->SegList)[k].nsyn;
        if ( nsyn > 0 ) {
            seg = &(b->SegList)[k];
            seg->vlam = (double *) malloc( (nsyn+1)*sizeof(double) );
            seg->diag = (double *) malloc( nsyn*sizeof(double) );
            seg->hseg = b->hseg;
            fac = (seg->hseg)/(pi*(seg->rp)*(seg->rd)*GA);
            vnew = 0.0;
            syn = seg->SynList;
            for ( j=0 ; j<nsyn ; j++ ) {
                (seg->vlam)[j] = vnew;
                vold = vnew;
                vnew = syn->vlam;
                (seg->diag)[j] = fac*(vnew-vold);
                syn = syn->NextSyn;
            }
            seg->vlam[nsyn] = vnew;

/*  Phase 2. - Allocate phi1, phi2, phi3 and diag */
            seg->phi1 = (double *) malloc( nsyn*sizeof(double) );
            seg->phi2 = (double *) malloc( nsyn*sizeof(double) );
            seg->phi3 = (double *) malloc( nsyn*sizeof(double) );
            seg->soln = (double *) malloc( (nsyn+1)*sizeof(double) );
            seg->vec = (double *) malloc( (nsyn+1)*sizeof(double) );
            seg->phi = (double *) malloc( (nsyn+1)*sizeof(double) );
            syn = seg->SynList;
            for ( j=0 ; j<nsyn ; j++ ) {
                (seg->phi1)[j] = 1.0-(syn->vlam);
                (seg->phi2)[j] = syn->vlam;
                (seg->phi3)[j] = -(syn->vsyn);
                syn = syn->NextSyn;
            }

/*  Phase 3. - Allocate and initialise conductances and firing times */
            seg->gold = (double *) malloc( nsyn*sizeof(double) );
            seg->gnew = (double *) malloc( nsyn*sizeof(double) );
            syn = seg->SynList;
            j = 0;
            while ( syn ) {
                (seg->gold)[j] = (seg->gnew)[j] = 0.0;
                interval = -(((double) NT*1000)/RATE)*log(ran(&ix, &iy, &iz));
                if ( fmod(interval,1.0) <= 0.5 ) {
                    syn->SynTime = -((int) floor(interval));
                } else {
                    syn->SynTime = -((int) ceil(interval));
                }
                syn = syn->NextSyn;
                j++;
            }

/*  Phase 4. - Initialise proximal and distal currents */
            seg->NewProx1 = 0.0;
            seg->NewProx2 = 0.0;
            seg->NewProx3 = 0.0;
            seg->NewDist1 = 0.0;
            seg->NewDist2 = 0.0;
            seg->NewDist3 = 0.0;
        }
    }
    return;
}


 /***********************************************
       Function to update status of synapses
  ***********************************************/
void Update_Synapses( branch *b )
{
    int i, j, k, nsyn;
    double interval, tmp;
    synapse *syn;
    segment *seg;

    if ( b->child ) Update_Synapses( b->child );
    if ( b->peer ) Update_Synapses( b->peer );

    for ( k=0 ; k<b->nseg ; k++ ) {
        seg = &(b->SegList)[k];

/*  Phase 1. - Update synaptic conductances */
        syn = seg->SynList;
        j = 0;
        while ( syn ) {
            (seg->gold)[j] = (seg->gnew)[j];
            (syn->SynTime)++;
            if ( syn->SynTime < 1 ) {
                (seg->gnew)[j] = 0.0;
            } else if ( syn->SynTime < syn->MaxStep ) {
                (seg->gnew)[j] = (syn->SynCond->g)[syn->SynTime];
            } else {
                (seg->gnew)[j] = 0.0;
                interval = -(((double) NT*1000)/RATE)*log(ran(&ix, &iy, &iz));
                if ( fmod(interval,1.0) <= 0.5 ) {
                    syn->SynTime = -((int) floor(interval));
                } else {
                    syn->SynTime = -((int) ceil(interval));
                }
            }
            syn = syn->NextSyn;
            j++;
        }

/*  Phase 2. - Update currents at segment endpoints */
        nsyn = seg->nsyn;
        if ( nsyn > 0 ) {
            seg->OldProx1 = seg->NewProx1;
            seg->OldProx2 = seg->NewProx2;
            seg->OldProx3 = seg->NewProx3;
            seg->OldDist1 = seg->NewDist1;
            seg->OldDist2 = seg->NewDist2;
            seg->OldDist3 = seg->NewDist3;

/*  Phase 3a. - Coefficient of VP */
            for ( j=0 ; j<nsyn ; j++ ) {
                (seg->phi)[j] = (seg->gnew)[j]*(seg->phi1)[j];
            }
            (seg->phi)[nsyn] = 0.0;
            solve( nsyn, seg->vlam, seg->phi, seg->soln);

/*  Phase 3b. - Repeat once */
            for ( j=0 ; j<nsyn ; j++ ) {
                for ( tmp=0.0,i=0 ; i<=j ; i++ ) {
                    tmp += (seg->diag)[i]*(seg->soln)[i];
                }
                (seg->phi)[j] -= tmp*(seg->gnew)[j];
            }
            solve( nsyn, seg->vlam, seg->phi, seg->soln);
            seg->NewProx1 = (seg->soln)[0];
            seg->NewDist1 = (seg->soln)[nsyn];

/*  Phase 4a. - Coefficient of VD */
            for ( j=0 ; j<nsyn ; j++ ) {
                (seg->phi)[j] = (seg->gnew)[j]*(seg->phi2)[j];
            }
            (seg->phi)[nsyn] = 0.0;
            solve( nsyn, seg->vlam, seg->phi, seg->soln);

/*  Phase 4b. - Repeat once */
            for ( j=0 ; j<nsyn ; j++ ) {
                for ( tmp=0.0,i=0 ; i<=j ; i++ ) {
                    tmp += (seg->diag)[i]*(seg->soln)[i];
                }
                (seg->phi)[j] -= tmp*(seg->gnew)[j];
            }
            solve( nsyn, seg->vlam, seg->phi, seg->soln);
            seg->NewProx2 = (seg->soln)[0];
            seg->NewDist2 = (seg->soln)[nsyn];

/*  Phase 5a. - Constant term */
            for ( j=0 ; j<nsyn ; j++ ) {
                (seg->phi)[j] = (seg->gnew)[j]*(seg->phi3)[j];
            }
            (seg->phi)[nsyn] = 0.0;
            solve( nsyn, seg->vlam, seg->phi, seg->soln);

/*  Phase 5b. - Repeat once */
            for ( j=0 ; j<nsyn ; j++ ) {
                for ( tmp=0.0,i=0 ; i<=j ; i++ ) {
                    tmp += (seg->diag)[i]*(seg->soln)[i];
                }
                (seg->phi)[j] -= tmp*(seg->gnew)[j];
            }
            solve( nsyn, seg->vlam, seg->phi, seg->soln);
            seg->NewProx3 = (seg->soln)[0];
            seg->NewDist3 = (seg->soln)[nsyn];
        }
    }
    return;
}


 /***************************************************
        Function to constuct sparse matrices
  ***************************************************/
void Update_Sparse_Matrices( branch *b, int *counter)
{
    int k, CurrentNode, nc;
    extern double pi, dt, *NewAmp;
    extern SparseMatrix lhs, rhs;
    branch *btmp;
    segment NewSeg, OldSeg;
    double SumL, SumR, tmp;

/* Step 1 - Recurse to the end of the dendrite */
    if ( b->child ) Update_Sparse_Matrices( b->child, counter);
    if ( b->peer ) Update_Sparse_Matrices( b->peer, counter);

/* Step 2 - Build matrix entries for distal node of branch */
    nc = b->nseg;
    tmp = 0.5*dt;
    CurrentNode = (b->first)-(nc-1);
    if ( b->child ) {
        btmp = b->child;
        SumR = SumL = 0.0;
        while ( btmp ) {
            NewSeg = (btmp->SegList)[0];
            lhs.a[*counter] = tmp*(NewSeg.NewProx2);
            rhs.a[*counter] = -tmp*(NewSeg.OldProx2);
            NewAmp[CurrentNode] += NewSeg.NewProx3;
            SumL += NewSeg.NewProx1;
            SumR += NewSeg.OldProx1;
            (*counter)++;
            btmp = btmp->peer;
        }
        NewSeg = (b->SegList)[nc-1];
        lhs.a[*counter] = tmp*(SumL-NewSeg.NewDist2);
        rhs.a[*counter] = -tmp*(SumR-NewSeg.OldDist2);
        NewAmp[CurrentNode] -= NewSeg.NewDist3;
        (*counter)++;
        lhs.a[*counter] = -tmp*(NewSeg.NewDist1);
        rhs.a[*counter] = tmp*(NewSeg.OldDist1);
        (*counter)++;
    } else {
        NewSeg = (b->SegList)[nc-1];
        lhs.a[*counter] = -tmp*(NewSeg.NewDist2);
        rhs.a[*counter] = tmp*(NewSeg.OldDist2);
        NewAmp[CurrentNode] = -NewSeg.NewDist3;
        (*counter)++;
        lhs.a[*counter] = -tmp*(NewSeg.NewDist1);
        rhs.a[*counter] = tmp*(NewSeg.OldDist1);
        (*counter)++;
    }

/* Step 3 - Build matrix entries for internal nodes of branch */
    for ( k=nc-1 ; k>0 ; k-- ) {
        CurrentNode++;
        OldSeg = NewSeg;
        NewSeg = (b->SegList)[k-1];
        lhs.a[*counter] = tmp*(OldSeg.NewProx2);
        rhs.a[*counter] = -tmp*(OldSeg.OldProx2);
        (*counter)++;
        lhs.a[*counter] = tmp*(OldSeg.NewProx1-NewSeg.NewDist2);
        rhs.a[*counter] = -tmp*(OldSeg.OldProx1-NewSeg.OldDist2);
        NewAmp[CurrentNode] = OldSeg.NewProx3-NewSeg.NewDist3;
        (*counter)++;
        lhs.a[*counter] = -tmp*(NewSeg.NewDist1);
        rhs.a[*counter] = tmp*(NewSeg.OldDist1);
        (*counter)++;
    }
    return;
}


 /**********************************************************
     Multiplies sparse matrix a[ ][ ] with vector v[ ]
  **********************************************************/
void Matrix_Vector_Multiply( SparseMatrix *a, double *v , double *b)
{
    int i, j, k, n;

    n = a->n;
    for ( j=0 ; j<n ; j++) {
        k = a->StartRow[j+1];
        for( b[j]=0.0,i=(a->StartRow[j]) ; i<k ; i++ ) {
            b[j] += (a->a[i])*v[a->col[i]];
        }
    }
    return;
}


 /***********************************************
        Allocate memory to a sparse matrix
  ***********************************************/
void Matrix_Malloc( SparseMatrix *a, int n, int w)
{
    a->a = (double *) malloc( w*sizeof(double) );
    a->col = (int *) malloc( w*sizeof(int) );
    a->StartRow = (int *) malloc( (n+1)*sizeof(int) );
    a->n = n;
    return;
}


 /**********************************************
     De-allocates memory of a sparse matrix
  **********************************************/
void Matrix_Free( SparseMatrix *a)
{
    free(a->a);
    free(a->col);
    free(a->StartRow);
    free(a);
}


void solve( int nsyn, double *vlam, double *b, double *soln)
{
    double sum;
    int k;

/*  Step 1. - Forward substitution phase */
    for ( sum=0.0,k=0 ; k<nsyn ; k++ ) {
        soln[k] = b[k];
        sum += vlam[k+1]*soln[k];
    }
    soln[nsyn] = b[nsyn]-sum;

/*  Step 2. - Backward substitution phase */
    for ( k=nsyn-1 ; k>=0 ; k-- ) soln[k] += soln[k+1];
    return;
}



 /************************************************************
         Function returns primitive uniform random number.
  ************************************************************/
double ran(unsigned int *ix, unsigned int *iy, unsigned int *iz)
{
    double tmp;

/*  1st item of modular arithmetic  */
    *ix = (171*(*ix))%30269;
/*  2nd item of modular arithmetic  */
    *iy = (172*(*iy))%30307;
/*  3rd item of modular arithmetic  */
    *iz = (170*(*iz))%30323;
/*  Generate random number in (0,1) */
    tmp = ((double) (*ix))/30269.0+((double) (*iy))/30307.0
          +((double) (*iz))/30323.0;
    return fmod(tmp,1.0);
}


 /********************************************************************
                    ALPHA for ACTIVATION OF SODIUM
  *******************************************************************/
double alfa_m( double volt )
{
    double tmp;
    static double fac;
    static int start=1;

    if ( start ) {
        fac = pow(3.0,0.1*CELSIUS-0.63);
        start = !start;
    }
    tmp = -0.1*(volt-25.0);
    if ( fabs(tmp)<0.001 ) {
        tmp = 1.0/(((tmp/24.0+1.0/6.0)*tmp+0.5)*tmp+1.0);
    } else {
        tmp = tmp/(exp(tmp)-1.0);
    }
    return tmp*fac;
}


 /********************************************************************
                    BETA for ACTIVATION OF SODIUM
  ********************************************************************/
double beta_m( double volt )
{
    double tmp;
    static double fac;
    static int start=1;

    if ( start ) {
        fac = pow(3.0,0.1*CELSIUS-0.63);
        start = !start;
    }
    tmp = volt/18.0;
    return 4.0*fac*exp(-tmp);
}


 /********************************************************************
                    ALPHA for INACTIVATION OF SODIUM
  ********************************************************************/
double alfa_h( double volt )
{
    double tmp;
    static double fac;
    static int start=1;

    if ( start ) {
        fac = pow(3.0,0.1*CELSIUS-0.63);
        start = !start;
    }
    tmp = 0.05*volt;
    return 0.07*fac*exp(-tmp);
}


 /********************************************************************
                    BETA for INACTIVATION OF SODIUM
  ********************************************************************/
double beta_h( double volt )
{
    double tmp;
    static double fac;
    static int start=1;

    if ( start ) {
        fac = pow(3.0,0.1*CELSIUS-0.63);
        start = !start;
    }
    tmp = -0.1*(volt-30.0);
    return fac/(exp(tmp)+1.0);
}


 /********************************************************************
                    ALPHA for ACTIVATION OF POTASSIUM
  ********************************************************************/
double alfa_n( double volt )
{
    double tmp;
    static double fac;
    static int start=1;

    if ( start ) {
        fac = pow(3.0,0.1*CELSIUS-0.63);
        start = !start;
    }
    tmp = -0.1*(volt-10.0);
    if ( fabs(tmp)<0.001 ) {
        tmp = 0.1/(((tmp/24.0+1.0/6.0)*tmp+0.5)*tmp+1.0);
    } else {
        tmp = 0.1*tmp/(exp(tmp)-1.0);
    }
    return tmp*fac;
}


 /********************************************************************
                    BETA for ACTIVATION OF POTASSIUM
  ********************************************************************/
double beta_n( double volt )
{
    double tmp;
    static double fac;
    static int start=1;

    if ( start ) {
        fac = pow(3.0,0.1*CELSIUS-0.63);
        start = !start;
    }
    tmp = 0.0125*volt;
    return 0.125*fac*exp(-tmp);
}


 /********************************************
        Computes a conductance profile
  ********************************************/
cond *ConductanceProfile( double dt,   /* Integration time step */
                          double tau,  /* Synaptic time constant */
                          double tol,  /* Determines off-threshold for synapse */
                          double gmax  /* Maximum conductance (mS) */ )
{
    int k;
    cond *out;
    double tmp, told, tnew;

    out = (cond *) malloc( sizeof(cond) );
    out->dt = dt;
    out->tau = tau;
    out->tol = tol;
    out->gmax = gmax;

/* Iterate to find duration of pulse */
    tmp = 1.0-log(tol);
    tnew = tmp;
    do {
        told = tnew;
        tnew = tmp+log(told);
    } while ( fabs(tnew-told)>5.e-11 );
    out->n = ((int) ceil(tau*tnew/dt));
    out->g = (double *) malloc( (out->n)*sizeof(double) );
    out->g[0] = 0.0;
    for ( k=1 ; k<(out->n) ; k++ ) {
        tmp = dt*((double) k)/tau;
        out->g[k] = gmax*tmp*exp(1.0-tmp);
    }
    return out;
}


void cgs(int *getmem, SparseMatrix *a, double *b, double *x, double tol)
{
    long int i, k, n, repeat;
    static int start=1;
    double rho_old, rho_new, alpha, beta, sigma, err;
    static double *p, *q, *r, *u, *v, *rb, *y;

/* Step 1 - Check memory status */
    n = a->n;
    if ( start ) {
        *getmem = 1;
        start = 0;
    }
    if ( *getmem ) {
        if ( p ) free(p);
        p = (double *) malloc( n*sizeof(double) );
        if ( q ) free(q);
        q = (double *) malloc( n*sizeof(double) );
        if ( r ) free(r);
        r = (double *) malloc( n*sizeof(double) );
        if ( u ) free(u);
        u = (double *) malloc( n*sizeof(double) );
        if ( v ) free(v);
        v = (double *) malloc( n*sizeof(double) );
        if ( rb ) free(rb);
        rb = (double *) malloc( n*sizeof(double) );
        if ( y ) free(y);
        y = (double *) malloc( n*sizeof(double) );
        *getmem = 0;
    }

/* Step 2 - Initialise residual, p[ ] and q[ ] */
    Matrix_Vector_Multiply( a, x, r);
    for ( rho_old=0.0,i=0 ; i<n ; i++ ) {
        r[i] = b[i]-r[i];
        rho_old += r[i]*r[i];
        rb[i] = r[i];
        p[i] = 0.0;
        q[i] = 0.0;
    }
    if ( rho_old<tol*((double) n) ) return;

/* The main loop */
    rho_old = 1.0;
    do {

/* Compute scale parameter for solution update */
        for ( rho_new=0.0,i=0 ; i<n ; i++ ) rho_new += r[i]*rb[i];
        beta = rho_new/rho_old;

/* Update u[ ] and p[ ] */
        for ( i=0 ; i<n ; i++ ) {
            u[i] = r[i]+beta*q[i];
            p[i] = u[i]+beta*(q[i]+beta*p[i]);
        }

/* Update v[ ] and compute sigma */
        Matrix_Vector_Multiply( a, p, v);
        for ( sigma=0.0,i=0 ; i<n ; i++ ) sigma += rb[i]*v[i];

/* Compute alpha and update q[ ], v[ ] and x[ ] */
        if ( sigma == 0.0 ) {
            printf(" Trouble ");
            for ( i=0 ; i<n ; i++ ) {
                printf("\n%20.16lf",v[i]);
                getchar( );
            }
        }
        alpha = rho_new/sigma;
        for ( i=0 ; i<n ; i++ ) {
            q[i] = u[i]-alpha*v[i];
            v[i] = alpha*(u[i]+q[i]);
            x[i] += v[i];
        }

 /* Update r[ ] and estimate error */
        Matrix_Vector_Multiply( a, v, y);
        for ( err=0.0,i=0 ; i<n ; i++ ) {
            r[i] -= y[i];
            err += r[i]*r[i];
        }
        rho_old = rho_new;
              if ( rho_old == 0.0 ) {
            printf(" Trouble rho_old ");
            for ( i=0 ; i<n ; i++ ) {
                printf("\n%20.16lf",v[i]);
                getchar( );
            }
        }
        repeat = ( err > tol*((double) n) );
    } while ( repeat );
    return;
}

Loading data, please wait...