#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
typedef struct SparseMatrix_t
{
double *a;
int *col;
int *StartRow;
int n;
struct SparseMatrix_t *l;
struct SparseMatrix_t *u;
} SparseMatrix;
typedef struct contact_t
{
int id; /* Identifies contact type */
double xp; /* Location of contact */
double pl; /* Measurement of physical length (micron) */
double amp; /* Strength of contact */
int xl; /* Left hand node */
int xr; /* Right hand node */
double frac; /* Fraction of input to left hand node */
struct contact_t *next; /* Address of next contact */
} contact;
typedef struct soma_t
{
/* Static biophysical properties of soma */
double cs; /* Somal membrane capacitance (mu F/cm^2) */
double ga; /* Intracellular conductance of soma (mS/cm) */
double gs; /* Membrane conductance of soma (mS/cm^2) */
/* Contact information */
contact *conlist; /* List of contacts */
} soma;
typedef struct branch_t
{
/* Connectivity of branch */
struct branch_t *parent; /* Pointer to parent branch */
struct branch_t *child; /* Pointer to child branch */
struct branch_t *peer; /* Pointer to a peer branch */
/* Physical properties of branch */
int nd; /* Number of nodes */
double xs; /* X-coord of somal endpoint */
double ys; /* Y-coord of somal endpoint */
double zs; /* Z-coord of somal endpoint */
double ds; /* Diameter of somal endpoint */
double xd; /* X-coord of distal endpoint */
double yd; /* Y-coord of distal endpoint */
double zd; /* Z-coord of distal endpoint */
double dd; /* Diameter of distal endpoint */
double *d; /* Diameter information */
double *x; /* Location information */
/* Biophysical properties of branch */
double cm; /* Dendritic membrane capacitance (mu F/cm^2) */
double ga; /* Intracellular conductance of dendrite (mS/cm) */
double gm; /* Membrane conductance of dendrite (mS/cm^2) */
/* 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 */
contact *conlist; /* List of contacts */
} branch;
typedef struct dendrite_t
{
branch *root; /* Pointer to root branch of dendrite */
double plen; /* Total length of dendrite */
} dendrite;
typedef struct neuron_t
{
int ndend; /* Number of dendrites */
dendrite *dendlist; /* Pointer to an array of dendrites */
soma *s; /* Soma structure */
} neuron;
/* Function type declarations */
neuron *LoadTestNeuron(char *);
void BuildTestDendrite( branch **, branch *),
RemoveBranch( branch **, branch *),
DestroyTestNeuron( neuron * ),
DestroyTestDendrite( branch *);
void BuildContactInfo(contact *, branch *, branch **);
int CountBranches( branch *, branch *),
CountContacts( branch *, branch *);
double branch_length( branch *, branch *);
int count_terminal_branches( branch *, branch *);
void output_properties( branch * );
void enumerate_nodes( branch *, int *);
void FreeSparseMatrix( SparseMatrix *),
MatrixVectorMultiplication( SparseMatrix *, double *, double *),
SparseMatrixMalloc( SparseMatrix *, int, int),
LU_Factor( SparseMatrix *, int *),
LU_Solve( SparseMatrix *, double *, double *);
void Generate_Dendrite(branch *, int *);
void Input_Current( branch *);
void Assign_Current( branch *, double *, double );
void count_dendritic_length ( branch *, double *);
void phys_lengths( branch *, double * );
/* Declaration of HH coefficient functions */
double alfa_h( double );
double alfa_m( double );
double alfa_n( double );
double beta_h( double );
double beta_m( double );
double beta_n( double );
/* Global definitions */
#define CELSIUS 18.5 /* Celsius temperature of neuron */
#define CS 1.0
#define GA 14.286
#define CM 1.0
#define GM 0.091
#define OUTPUT "contactinfo.dat"
#define TEND 2.0
#define NT 100
#define DT 0.1
#define SIN 0.0e-3
#define T 100.0
#define RS 0.002
#define NNODES 100
#define SENTRY 50.0
/* Global Variables */
SparseMatrix lhs, rhs, lhs1;
int main( int argc, char **argv )
{
int i, k, id, start, nodes, nc, in, FirstNode, counter,
nb, nspk, first;
double jo, je, jn, ho, he, hn, mo, me, mn, no, ne, nn,
vo, ve, vn, xo, xn;
double pi, as;
double *v, *x, max, gama, vs, frac, arg,
sum, tmp, tmp1, tmp2, dt, tnow, tout, left,
rite, gval, cval, len, h, vmid;
double g_na = 120.0, g_k = 36.0, g_l = 0.3, v_na = 55.0,
v_k = -72.0, v_l = -49.416, fact=25.0;
neuron *n;
extern SparseMatrix lhs, rhs, lhs1, ptr;
branch *bnow;
FILE *fp;
/* Compute ancillary variables */
pi = 4.0*atan(1.0);
as = 4.0*pi*RS*RS;
/* Compute Equilibrium Potential and equilibrium states */
vo = -62.0;
mo = alfa_m(vo)/(alfa_m(vo)+beta_m(vo));
ho = alfa_h(vo)/(alfa_h(vo)+beta_h(vo));
no = alfa_n(vo)/(alfa_n(vo)+beta_n(vo));
jo = g_na*pow(mo,3)*ho*(vo-v_na)+g_k*pow(no,4)*(vo-v_k)+g_l*(vo-v_l);
vn = -58.0;
mn = alfa_m(vn)/(alfa_m(vn)+beta_m(vn));
hn = alfa_h(vn)/(alfa_h(vn)+beta_h(vn));
nn = alfa_n(vn)/(alfa_n(vn)+beta_n(vn));
jn = g_na*pow(mn,3)*hn*(vn-v_na)+g_k*pow(nn,4)*(vn-v_k)+g_l*(vn-v_l);
if ( jo*jn > 0.0 ) {
printf(" No zero found \n");
return(0);
} else {
while ( fabs(vo-vn) > 5.e-7 ) {
ve = 0.5*(vo+vn);
me = alfa_m(ve)/(alfa_m(ve)+beta_m(ve));
he = alfa_h(ve)/(alfa_h(ve)+beta_h(ve));
ne = alfa_n(ve)/(alfa_n(ve)+beta_n(ve));
je = g_na*pow(me,3)*heq*(ve-v_na)+g_k*pow(ne,4)*(ve-v_k)+g_l*(veq-v_l);
if ( je*jo > 0.0 ) {
vo = ve;
} else {
vn = ve;
}
}
}
vn = ve;
mn = alfa_m(vn)/(alfa_m(vn)+beta_m(vn));
hn = alfa_h(vn)/(alfa_h(vn)+beta_h(vn));
nn = alfa_n(vn)/(alfa_n(vn)+beta_n(vn));
v_na -= ve;
v_k -= ve;
v_l -= ve;
/* Load sampled neuron */
if ( argc != 2 ) {
printf("\n Invoke program with load <input>\n");
return(1);
} else {
n = LoadTestNeuron( argv[1] );
if ( !n ) {
printf("\n Failed to find test neuron\n");
return(1);
}
len = 0.0;
nb = 0;
for ( i=0 ; i<n->ndend ; i++) count_branch( n->dendlist[i].root, &nb );
for ( i=0 ; i<n->ndend ; i++) count_dendritic_length( n->dendlist[i].root, &length );
h = len/((double) NNODES - nb);
for( i=0 ; i<n->ndend ; i++) phys_lengths( n->dendlist[i].root, &h );
}
/* Enumerate Nodes */
FirstNode = 0;
for ( k=0 ; k<n->ndend ; k++ ) enumerate_nodes( n->dendlist[k].root, &FirstNode );
for ( k=0 ; k<n->ndend ; k++ ) n->dendlist[k].root->junct = FirstNode;
printf("Number of nodes is %d\n", FirstNode+1);
/* Construct Sparse Matrices */
counter = 0;
nodes = FirstNode+1;
mat_malloc( &lhs, nodes, 3*nodes-2 );
mat_malloc( &rhs, nodes, 3*nodes-2 );
mat_malloc( &lhs1, nodes, 3*nodes-2);
lhs.n = rhs.n = nodes;
lhs.a[3*nodes-3] = rhs.a[3*nodes-3] = 0.0;
lhs.start_row[0] = rhs.start_row[0] = 0;
for ( k=0 ; k<n->ndend ; k++ ) {
bnow = n->dendlist[k].root;
Generate_Dendrite( bnow, &counter);
}
for ( k=0 ; k<n->ndend ; k++ ) {
bnow = n->dendlist[k].root;
lhs.a[counter] = 0.5*bnow->d[1]*bnow->pl[1];
rhs.a[counter] = -(bnow->d[0])*(bnow->d[1])/bnow->pl[1];
lhs.col[counter] = rhs.col[counter] = bnow->first;
lhs.a[3*nodes-3] += 1.5*(bnow->d[0])*bnow->pl[1];
rhs.a[3*nodes-3] += (bnow->d[0])*(bnow->d[1])/bnow->pl[1];
counter++;
}
lhs.col[3*nodes-3] = rhs.col[3*nodes-3] = nodes-1;
lhs.start_row[nodes] = rhs.start_row[nodes] = 3*nodes-2;
for( i=0; i<n->ndend ; i++ ) Input_Current(n->dendlist[i].root);
dt = 1.0/((double) NT);
for ( i=0 ; i<3*nodes-2 ; i++ ) {
rhs.a[i] *= GA;
rhs.a[i] += GM*lhs.a[i];
lhs.a[i] *= CM;
rhs.a[i] *= 0.5*dt;
lhs.a[i] += rhs.a[i];
rhs.a[i] = lhs.a[i] - 2.0*rhs.a[i];
}
left = lhs.a[3*nodes-3]+(4.0*as/pi)*CS;
rite = rhs.a[3*nodes-3]+(4.0*as/pi)*CS;
v = (double *) malloc( (nodes)*sizeof(double) );
x = (double *) malloc( (nodes)*sizeof(double) );
for ( i=0 ; i<nodes ; i++ ) v[i] = 0.0;
/* Initialise temporal integration */
tnow = 0.0;
tout = DT;
vold = vnew = vmid = 0.0;
nspk = 0;
first = 1;
while ( tnow < TEND ) {
tnow += dt;
vs = veq+v[nodes-1];
mold = mnew;
nold = nnew;
hold = hnew;
tmp1 = dt*alfa_m(vs);
tmp2 = dt*beta_m(vs);
tmp = tmp1+tmp2;
mnew = mold+(tmp1-mold*tmp)/(1.0+0.5*tmp);
tmp1 = dt*alfa_n(vs);
tmp2 = dt*beta_n(vs);
tmp = tmp1+tmp2;
nnew = nold+(tmp1-nold*tmp)/(1.0+0.5*tmp);
tmp1 = dt*alfa_h(vs);
tmp2 = dt*beta_h(vs);
tmp = tmp1+tmp2;
hnew = hold+(tmp1-hold*tmp)/(1.0+0.5*tmp);
gval = g_na*hnew*pow(mnew,3);
cval = v_na*gval;
tmp = g_k*pow(nnew,4);
gval += tmp;
cval += tmp*v_k;
gval += g_l;
cval += g_l*v_l;
tmp = as*4.0*dt/pi;
gval *= 0.5*tmp*fact;
cval *= tmp*fact;
lhs.a[3*nodes-3] = left + gval;
rhs.a[3*nodes-3] = rite - gval;
for ( k=0 ; k<3*nodes-2 ; k++ ) {
lhs1.a[k] = lhs.a[k];
lhs1.col[k] = lhs.col[k];
}
for ( k=0 ; k<=nodes ; k++ ) lhs1.start_row[k] = lhs.start_row[k];
OldLU_Factor(&lhs);
fp = fopen("out-l2.dat","w");
for ( i=0 ; i<nodes ; i++ ) {
for ( k=lhs.l->start_row[i] ; k<lhs.l->start_row[i+1] ; k++ ) {
fprintf(fp,"%d \t (%d,%d) \t %lf\n", k, i, lhs.l->col[k], lhs.l->a[k]);
}
}
fclose(fp);
fp = fopen("out-u2.dat","w");
for ( i=0 ; i<nodes ; i++ ) {
for ( k=lhs.u->start_row[i] ; k<lhs.u->start_row[i+1] ; k++ ) {
fprintf(fp,"%d \t (%d,%d) \t %lf\n", k, i, lhs.u->col[k], lhs.u->a[k]);
}
}
fclose(fp);
//return 0;
LU_Factor(&lhs1, &first);
fp = fopen("out-l1.dat","w");
for ( i=0 ; i<nodes ; i++ ) {
for ( k=lhs1.l->start_row[i] ; k<lhs1.l->start_row[i+1] ; k++ ) {
fprintf(fp,"%d \t (%d,%d) \t %lf\n", k, i, lhs1.l->col[k], lhs1.l->a[k]);
}
}
fclose(fp);
fp = fopen("out-u1.dat","w");
for ( i=0 ; i<nodes ; i++ ) {
for ( k=lhs1.u->start_row[i] ; k<lhs1.u->start_row[i+1] ; k++ ) {
fprintf(fp,"%d \t (%d,%d) \t %lf\n", k, i, lhs1.u->col[k], lhs1.u->a[k]);
}
}
fclose(fp);
tmp = 0.0;
for ( i=0 ; i<nodes ; i++ ) {
for ( k=lhs.l->start_row[i] ; k<lhs.l->start_row[i+1] ; k++ ) {
if ( fabs(lhs1.l->a[k]-lhs.l->a[k]) > tmp ) {
tmp = fabs(lhs1.l->a[k]-lhs.l->a[k]);
}
}
}
if ( tmp > 5.e-12 ) {
printf("\n trouble %lf %lf",tnow, tmp);
getchar( );
}
// return 0;
mat_vec_mult(&rhs,v,x);
x[nodes-1] -= 4.0*dt*SIN/pi;
x[nodes-1] += cval;
tmp = 4.0*dt/pi;
if (tnow<T) for (i=0;i<n->ndend;i++) Assign_Current(n->dendlist[i].root,x,tmp);
LU_Solve( &lhs, v, x );
vnew = v[nodes-1];
if ( vmid > SENTRY ) {
if ( vold < vmid && vnew < vmid ) {
tmp1 = tnow + 0.5*dt*(vnew-vold)/(2.0*vmid-vold-vnew);
tmp2 = vmid + 0.125*pow(vnew-vold,2)/(2.0*vmid-vold-vnew);
if ( nspk == 0 ) fp = fopen("out.res", "w");
if ( nspk != 0 ) fp = fopen("out.res", "a");
nspk++;
fprintf(fp,"Spike %d of %lf mv at time %lf ms \n", nspk, tmp2, tmp1);
fclose(fp);
}
}
vold = vmid;
vmid = vnew;
if ( tnow > tout ) {
printf("\rReached time %5.1lf ms\t", tout);
printf("\nNumerical Voltage %12.6lf mV\n",v[nodes-1]);
tout += DT;
}
}
/* Count contacts
for ( n=k=0 ; k<n->ndend ; k++ ) {
nc += count_contacts( n->dendlist[k].root, n->dendlist[k].root);
}
printf("\n Located %d contacts on dendrites", nc);
printf("\n Located %d contacts on soma", n->s->ncon);
printf("\n"); */
DestroyTestNeuron( n );
return(0);
}
/*************************************************
Function To Load A Test Neuron
*************************************************/
neuron *LoadTestNeuron(char *filename)
{
int j, k, ncon, n, id, connected, ignored;
double tmp, piby2, xl, xr, dl, dr, px, py, pz, min, radius, dx;
neuron *cell;
contact *newcon;
branch *bold, *bnew, *FirstBranch;
char temp[100];
FILE *fp;
/* STEP 1. - Open neuron data file */
printf("\nOpening file %s\n",filename);
if ( (fp=fopen(filename,"r"))==NULL ) return NULL;
/* STEP 2. - Get memory for neuron structure */
cell = (neuron *) malloc( sizeof(neuron) );
/* STEP 3. - Get branch and contact data */
bo = NULL;
while ( fscanf(input,"%s", temp)!=EOF ) {
if ( strcmp(temp,"Branch") == 0 || strcmp(temp,"branch") == 0 ) {
fscanf(fp, "%d", &n);
printf("Found a branch\n");
bnew = (branch *) malloc( sizeof(branch) );
bnew->nd = n;
if ( bold ) {
bold->child = bnew;
} else {
FirstBranch = bnew;
}
bnew->parent = bold;
bnew->peer = NULL;
bnew->child = NULL;
/* STEP 3b. - Initialise branch */
bnew->d = (double *) malloc( n*sizeof(double) );
bnew->x = (double *) malloc( n*sizeof(double) );
bnew->cm = CM;
bnew->gm = GM;
bnew->ga = GA;
bnew->conlist = NULL;
/* STEP 3c. - Read branch morphology */
fscanf(fp,"%lf %lf %lf %lf", &(bnew->xs), &(bnew->ys), &(bnew->zs), &(bnew->ds) );
fscanf(fp,"%lf %lf %lf %lf", &(bnew->xd), &(bnew->yd), &(bnew->zd), &(bnew->dd) );
fscanf(fp,"%lf", &len );
dx = len/((double) n-1 );
dd = (bnew->dd-bnew->ds)/((double) n-1);
for ( j=0 ; j<n ; j++ ) {
bnew->x[j] = dx*((double) j);
bnew->d[j] = ds+dd*((double) j);
}
bold = bnew;
} else if ( strcmp(temp, "Marker") == 0 || strcmp(temp, "marker") == 0 ) {
/* STEP 3d. - Initialise marker */
printf("Found and initialised a branch contact\n");
newcon = (contact *) malloc( sizeof(contact) );
newcon->next = NULL;
fscanf(fp,"%lf %lf", &newcon->xp, &newcon->amp );
if ( bnew->conlist == NULL ) {
bnew->conlist = newcon;
} else {
oldcon = bnew->conlist;
while ( oldcon->next ) oldcon = oldcon->next;
oldcon->next = newcon;
}
} else {
printf("Unrecognised dendritic feature\n");
}
}
fclose(fp);
/* STEP 4. - Count dendritic branches at soma */
bold = FirstBranch;
n = 0;
while ( bold ) {
bnew = FirstBranch;
do {
tmp = pow(bold->xs-bnew->xd,2)+pow(bold->ys-bnew->yd,2)
+pow(bold->zs-bnew->zd,2);
connected = ( tmp < 0.01 );
bnew = bnew->child;
} while ( bnew && !connected );
if ( !connected ) n++;
bold = bold->child;
}
cell->ndend = n;
printf("\n\nTree contains %d individual dendrite(s) ...\n", n);
/* STEP 5. - Identify somal dendrites but extract nothing */
cell->dendlist = (dendrite *) malloc( (cell->ndend)*sizeof(dendrite) );
bold = FirstBranch;
n = 0;
while ( n < cell->ndend ) {
bnew = FirstBranch;
do {
tmp = pow(bold->xs-bnew->xd,2)+pow(bold->ys-bnew->yd,2)
+pow(bold->zs-bnew->zd,2);
connected = ( tmp < 0.01 );
bnew = bnew->child;
} while ( bnew && !connected );
if ( !connected ) cell->dendlist[n++].root = bold;
bold = bold->child;
}
/* STEP 6. - Extract root of each dendrite from dendrite list */
for ( k=0 ; k<cell->ndend ; k++ ) {
bold = cell->dendlist[k].root;
RemoveBranch( &FirstBranch, bold);
}
/* STEP 7. - Build each test dendrite from its root branch */
for ( k=0 ; k<cell->ndend ; k++ ) {
BuildTestDendrite( &FirstBranch, cell->dendlist[k].root);
}
if ( FirstBranch != NULL ) printf("\nWarning: Unconnected branch segments still exist\n");
return cell;
}
/**************************************************
Function to remove a branch from a branch list
**************************************************/
void RemoveBranch( 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 build a test dendrite from its root
********************************************************/
void BuildTestDendrite( 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;
/* Search if proximal end of bnow is connected to distal end of root */
tmp = pow(bnow->xs-root->xd,2)+pow(bnow->ys-root->yd,2)+
pow(bnow->zs-root->zd,2);
if ( tmp <= 0.01 ) {
/* Take bnow out of 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 ) BuildTestDendrite( head, root->child);
if ( root->peer ) BuildTestDendrite( head, root->peer);
return;
}
/*****************************************************
Function to destroy a NEURON
*****************************************************/
void DestroyTestNeuron( neuron *cell)
{
int k;
contact *prevcon, *nextcon;
/* Free Soma */
if ( cell->s != NULL ) {
prevcon = cell->s->conlist;
while ( prevcon ) {
nextcon = prevcon->next;
free ( prevcon );
prevcon = nextcon;
}
free ( cell->s );
}
for ( k=0 ; k<cell->ndend ; k++ ) DestroyTestDendrite( cell->dendlist[i].root );
free(cell);
return;
}
/************************************************
Function to destroy Test DENDRITE
************************************************/
void DestroyTestDendrite( branch *b )
{
contact *prevcon, *nextcon;
if ( b->child ) DestroyTestDendrite(b->child);
if ( b->peer ) DestroyTestDendrite(b->peer);
free( b->x );
free( b->d );
prevcon = b->conlist;
while ( prevcon ) {
nextcon = prevcon->next;
free ( prevcon );
prevcon = nextcon;
}
free( b );
return;
}
/****************************************************
Function to count number of branches
from current branch to dendritic tip
****************************************************/
int CountBranches( branch *bstart, branch *bnow)
{
static int n;
if ( bstart == bnow ) n = 0;
if ( bnow ) {
if ( bnow->child ) CountBranches( bstart, bnow->child);
if ( bnow->peer ) CountBranches( bstart, bnow->peer);
n++;
}
return n;
}
/*****************************************************
Function to count number of contacts
from current branch to the dendritic tip.
*****************************************************/
int CountContacts( branch *bstart, branch *bnow)
{
static int n;
contact *con;
if ( bstart == bnow ) n = 0;
if ( bnow ) {
if ( bnow->child ) CountContacts(bstart, bnow->child);
if ( bnow->peer ) CountContacts(bstart, bnow->peer);
con = bnow->conlist;
while ( con ) {
n++;
con = con->next;
}
}
return n;
}
/******************************************************
Function to constuct sparse matrices
******************************************************/
void Generate_Dendrite( branch *bnow, int *counter)
{
int i, k;
extern sparse_mat lhs, rhs;
branch *btmp;
double SumL, SumR;
/* Step 1 - Recurse to the end of the dendrite */
if ( bnow->child != NULL ) Generate_Dendrite( bnow->child, counter);
if ( bnow->peer != NULL ) Generate_Dendrite( bnow->peer, counter);
for ( k=bnow->first-bnow->nobs+2,i=bnow->nobs-1 ; i>0 ; i--,k++ ) {
/* Step 2 - Fill in matrix entries for terminal points */
if ( bnow->child == NULL && i == bnow->nobs - 1 ) {
lhs.a[*counter] = 1.5*(bnow->d[i])*(bnow->pl[i]-bnow->pl[i-1]);
rhs.a[*counter] = (bnow->d[i-1]*bnow->d[i])/(bnow->pl[i]-bnow->pl[i-1]);
lhs.col[*counter] = rhs.col[*counter] = k;
(*counter)++;
lhs.a[*counter] = 0.5*(bnow->d[i-1])*(bnow->pl[i] - bnow->pl[i-1]);
rhs.a[*counter] = -(bnow->d[i-1]*bnow->d[i])/(bnow->pl[i] - bnow->pl[i-1]);
if ( k == bnow->first ) {
lhs.col[*counter] = rhs.col[*counter] = bnow->junct;
} else {
lhs.col[*counter] = rhs.col[*counter] = k + 1;
}
(*counter)++;
lhs.start_row[k+1] = rhs.start_row[k+1] = *counter;
/* Step 3 - Fill in matrix entries for branch points */
} else if ( bnow->child != NULL && i == bnow->nobs - 1 ) {
btmp = bnow->child;
SumR = SumL = 0.0;
while ( btmp != NULL ) {
lhs.a[*counter] = 0.5*btmp->d[1]*btmp->pl[1];
rhs.a[*counter] = -(btmp->d[0]*btmp->d[1])/btmp->pl[1];
lhs.col[*counter] = rhs.col[*counter] = btmp->first;
(*counter)++;
SumL += 1.5*btmp->d[0]*btmp->pl[1];
SumR += (btmp->d[0]*btmp->d[1])/btmp->pl[1];
btmp = btmp->peer;
}
lhs.a[*counter] = SumL+1.5*bnow->d[i]*(bnow->pl[i]-bnow->pl[i-1]);
rhs.a[*counter] = SumR+(bnow->d[i-1]*bnow->d[i])/(bnow->pl[i]-bnow->pl[i-1]);
lhs.col[*counter] = rhs.col[*counter] = k;
(*counter)++;
lhs.a[*counter] = 0.5*(bnow->d[i-1])*(bnow->pl[i]-bnow->pl[i-1]);
rhs.a[*counter] = -(bnow->d[i-1]*bnow->d[i])/(bnow->pl[i] - bnow->pl[i-1]);
if ( k == bnow->first ) {
lhs.col[*counter] = rhs.col[*counter] = bnow->junct;
} else {
lhs.col[*counter] = rhs.col[*counter] = k + 1;
}
(*counter)++;
lhs.start_row[k+1] = rhs.start_row[k+1] = *counter;
} else {
/* Step 4 - Fill in matrix entries for internal point */
lhs.a[*counter] = 0.5*(bnow->d[i+1])*(bnow->pl[i+1] - bnow->pl[i]);
rhs.a[*counter] = -(bnow->d[i]*bnow->d[i+1])/(bnow->pl[i+1] - bnow->pl[i]);
lhs.col[*counter] = rhs.col[*counter] = k - 1;
(*counter)++;
lhs.a[*counter] = 1.5*(bnow->d[i])*(bnow->pl[i+1] - bnow->pl[i-1]);
rhs.a[*counter] = (bnow->d[i-1]*bnow->d[i])/(bnow->pl[i] - bnow->pl[i-1])
+ (bnow->d[i]*bnow->d[i+1])/(bnow->pl[i+1] - bnow->pl[i]);
lhs.col[*counter] = rhs.col[*counter] = k;
(*counter)++;
lhs.a[*counter] = 0.5*(bnow->d[i-1])*(bnow->pl[i] - bnow->pl[i-1]);
rhs.a[*counter] = -(bnow->d[i-1]*bnow->d[i])/(bnow->pl[i] - bnow->pl[i-1]);
lhs.col[*counter] = rhs.col[*counter] = k + 1;
if ( k == bnow->first ) {
lhs.col[*counter] = rhs.col[*counter] = bnow->junct;
} else {
lhs.col[*counter] = rhs.col[*counter] = k + 1;
}
(*counter)++;
lhs.start_row[k+1] = rhs.start_row[k+1] = *counter;
}
}
return;
}
/***************************************************************
Function to build CONTACT information
***************************************************************/
void BuildContactInfo(contact *con, branch *b, branch **bopt)
{
int k;
double px, py, pz, tmp, xold, xnew, yold, ynew, zold, znew,
numer, denom, xmin, ymin, zmin, min;
px = con->xc;
py = con->yc;
pz = con->zc;
/* STEP 1. - First stage is different from others */
xnew = b->x[0]; ynew = b->y[0]; znew = b->z[0];
min = sqrt(pow(xnew-px,2)+pow(ynew-py,2)+pow(znew-pz,2));
if ( con->sd == NULL || ( con->sd != NULL && min < con->sd ) ) {
con->sd = min;
con->xp = xnew; con->yp = ynew; con->zp = znew;
con->pl = 0.0;
*bopt = b;
}
/* STEP 2. - Second stage compares points and projected points */
for ( k=1 ; k<b->nobs ; k++ ) {
xold = xnew; yold = ynew; zold = znew;
xnew = b->x[k]; ynew = b->y[k]; znew = b->z[k];
numer = (xnew-xold)*(px-xold)+(ynew-yold)*(py-yold)+(znew-zold)*(pz-zold);
denom = pow(xnew-xold,2)+pow(ynew-yold,2)+pow(znew-zold,2);
/* STEP 2a. - Project onto branch */
if ( 0.0 <= numer && numer <= denom ) {
tmp = numer/denom;
xmin = (1.0-tmp)*xold+tmp*xnew;
ymin = (1.0-tmp)*yold+tmp*ynew;
zmin = (1.0-tmp)*zold+tmp*znew;
min = sqrt(pow(xmin-px,2)+pow(ymin-py,2)+pow(zmin-pz,2));
if ( !(con->sd) || ( con->sd && min < con->sd ) ) {
con->sd = min;
con->xp = xmin; con->yp = ymin; con->zp = zmin;
con->pl = (1.0-tmp)*b->pl[k-1]+tmp*b->pl[k];
*bopt = b;
}
}
/* STEP 2b. - Check proximity to points of branch */
min = sqrt(pow(xnew-px,2)+pow(ynew-py,2)+pow(znew-pz,2));
if ( !(con->sd) || ( con->sd && min < con->sd ) ) {
con->sd = min;
con->xp = xnew; con->yp = ynew; con->zp = znew;
con->pl = b->pl[k];
*bopt = b;
}
}
return;
}
/***************************************************************
Function to find length of dendrite from
current branch to tips.
****************************************************************/
double branch_length( branch *bstart, branch *bnow)
{
static double length;
if ( bstart == bnow ) length = 0.0;
if ( bnow ) {
if ( bnow->child ) branch_length(bstart, bnow->child);
if ( bnow->peer ) branch_length(bstart, bnow->peer);
length += bnow->p_len;
}
return length;
}
/****************************************************************
Function to count number of terminal branches
****************************************************************/
int count_terminal_branches( branch *bstart, branch *bnow)
{
static int n;
if ( bstart == bnow ) n = 0;
if ( bnow ) {
if ( bnow->child ) count_terminal_branches(bstart, bnow->child);
if ( bnow->peer ) count_terminal_branches(bstart, bnow->peer);
if ( !bnow->child ) n++;
}
return n;
}
/****************************************************************
Function to output branch diameters
****************************************************************/
void output_properties( branch *b )
{
int i, k;
static int start=1;
double dold, dnew, len, xold, yold, zold, xnew, ynew, znew, dx, dy, dz, size;
branch *bran;
FILE *fp;
if ( b->child ) output_properties(b->child);
if ( b->peer ) output_properties(b->peer);
if ( start ) {
fp = fopen("output","w");
start = 0;
} else {
fp = fopen("output","a");
fprintf(fp,"\n");
}
/* Outputs branch lengths, diameters, surface areas etc.
for ( k=0 ; k<b->nobs ; k++ ) {
fprintf(fp,"%6.2lf \t %6.2lf \t %6.2lf \t %6.2lf \n", b->pl[k], b->d[k], b->sa[k], b->el[k]);
}*/
/* Decomposes branches into lengths of uniform diameter
len = xold = b->pl[1];
dold = b->d[1];
for ( k=2 ; k<b->nobs ; k++ ) {
xnew = b->pl[k];
dnew = b->d[k];
if ( dnew != dold ) {
len += 0.5*(xnew-xold);
fprintf(fp,"%6.2lf \t %6.2lf \n", len, dold);
len = 0.5*(xnew-xold);
} else {
len += xnew-xold;
}
xold = xnew;
dold = dnew;
}
fprintf(fp,"%6.2lf \t %6.2lf \n", len, dold); */
/* Constructs diameters of a branch and its children/peers */
if ( b->child ) {
fprintf(fp,"%6.2lf \t %6.2lf \t", b->d[(b->nobs)-1], b->child->d[1]);
bran = b->child;
while ( bran->peer ) {
bran = bran->peer;
fprintf(fp,"%6.2lf \t", bran->d[1]);
}
}
/* Prints out branch lengths
printf("\nBranch length %6.2lf, %6.2lf, %6.2lf", b->p_len, b->d[0], b->d[b->nobs-1] );
getchar( ); */
fclose(fp);
return;
}
/**********************************************************
Function to enumerate the nodes on a dendrite
**********************************************************/
void enumerate_nodes(branch *bnow, int *FirstNode )
{
branch *btmp;
if ( (bnow->child) != NULL ) enumerate_nodes( bnow->child, FirstNode );
if ( (bnow->peer) != NULL ) enumerate_nodes( bnow->peer, FirstNode );
if ( bnow->child != NULL ) {
btmp = bnow->child;
while( btmp != NULL ){
btmp->junct = *FirstNode;
btmp = btmp->peer;
}
}
bnow->first = *FirstNode + bnow->nobs - 2;
*FirstNode += (bnow->nobs)-1;
return;
}
/***************************************************
Allocate memory to a sparse matrix
***************************************************/
void SparseMatrixMalloc( 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;
a->l = malloc( sizeof(SparseMatrix) );
a->u = malloc( sizeof(SparseMatrix) );
a->l->a = (double *) malloc( (2*n-1)*sizeof(double) );
a->l->col = (int *) malloc( (2*n-1)*sizeof(int) );
a->l->StartRow = (int *) malloc( (n+1)*sizeof(int) );
a->l->n = n;
a->u->a = (double *) malloc( (2*n-1)*sizeof(double) );
a->u->col = (int *) malloc( (2*n-1)*sizeof(int) );
a->u->StartRow = (int *) malloc( (n+1)*sizeof(int) );
a->u->n = n;
return;
}
/********************************************************
Multiplies sparse matrix a[ ][ ] with vector v[ ]
********************************************************/
void MatrixVectorMultiplication( 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;
}
/**********************************************
De-allocates memory of a sparse matrix
**********************************************/
void FreeSparseMatrix( SparseMatrix *a)
{
free(a->a);
free(a->col);
free(a->StartRow);
free(a);
}
/***************************************************************
Function To Factorise A Sparse Matrix
***************************************************************/
void LU_Factor( SparseMatrix *m, int *start)
{
double tmp, sum;
int i, j, k, r, n, cl, cu, nrow;
/* Step 1. - Identify matrix dimension */
n = m->n;
/* Step 2. - Fill column vectors for triangular matrices */
if ( *start ) {
cl = cu = 0;
for ( i=k=0 ; i<n ; i++ ) {
m->l->StartRow[i] = cl;
m->u->StartRow[i] = cu;
while ( m->col[k] < i ) m->l->col[cl++] = m->col[k++];
m->l->col[cl++] = m->col[k];
m->u->col[cu++] = m->col[k++];
while ( k < m->StartRow[i+1] ) m->u->col[cu++] = m->col[k++];
}
m->l->StartRow[n] = cl;
m->u->StartRow[n] = cu;
*start = 0;
}
/* Step 3. - Fill first row of L and U */
m->l->a[0] = 1.0;
for ( k=0 ; k < m->u->StartRow[1] ; k++ ) m->u->a[k] = m->a[k];
/* Step 4. - Fill remaining entries row by row */
cl = 1;
k = cu = m->u->StartRow[1];
for ( i=1 ; i<n ; i++ ) {
while ( m->col[k] < i ) { // Fill lower matrix
sum = m->a[k];
for ( j=m->l->StartRow[i] ; j<cl ; j++ ) {
nrow = m->u->StartRow[m->l->col[j]];
while ( m->u->col[nrow] < m->col[k] ) nrow++;
sum -= (m->l->a[j])*(m->u->a[nrow]);
}
nrow = m->u->StartRow[m->l->col[cl]];
while ( m->u->col[nrow] < m->col[k] ) nrow++;
m->l->a[cl++] = sum/(m->u->a[nrow]);
k++;
}
m->l->a[cl++] = 1.0; // Diagonal entry of lower
while ( m->col[k] >= i && k < m->start_row[i+1] ) { // Fill upper matrix
sum = m->a[k];
for ( j=m->l->StartRow[i] ; j<m->l->StartRow[i+1]-1; j++ ) {
nrow = m->u->StartRow[m->l->col[j]];
while ( m->u->col[nrow] < m->col[k] && nrow < m->u->start_row[m->l->col[j]+1] ) nrow++;
sum -= (m->l->a[j])*(m->u->a[nrow]);
}
k++;
m->u->a[cu++] = sum;
}
}
return;
}
/****************************************************
Function to Solve the matrix problem
****************************************************/
void LU_Solve( SparseMatrix *m, double *x, double *b )
{
int i, j;
double *z;
z = (double *) malloc( (m->n)*sizeof(double) );
for ( i=0 ; i<m->n ; i++ ) {
z[i] = b[i];
for ( j=m->l->StartRow[i] ; j<m->l->StartRow[i+1]-1 ; j++ ) {
z[i] -= (m->l->a[j])*(z[(m->l->col[j])]);
}
z[i] /= m->l->a[m->l->StartRow[i+1]-1];
}
for ( i=(m->n)-1 ; i>=0 ; i-- ) {
x[i] = z[i];
for ( j=m->u->StartRow[i]+1 ; j<m->u->StartRow[i+1] ; j++ ) {
x[i] -= (m->u->a[j])*(x[m->u->col[j]]);
}
x[i] /= m->u->a[m->u->StartRow[i]];
}
free(z);
return;
}
/*************************************************************
Function to input current to dendrite
**************************************************************/
void Input_Current( branch *bnow )
{
int k;
double tmp;
if ( bnow->child != NULL ) Input_Current( bnow->child );
if ( bnow->peer != NULL ) Input_Current( bnow->peer );
if ( bnow->conlist != NULL ) {
if ( bnow->conlist->pl <= bnow->pl[1] ) {
bnow->conlist->xl = bnow->junct;
bnow->conlist->xr = bnow->first;
bnow->conlist->frac = (1.0-(bnow->conlist->pl/bnow->pl[1]));
} else {
k = 1;
while ( bnow->conlist->pl > bnow->pl[k] ) k++;
bnow->conlist->xl = bnow->first-k+2;
bnow->conlist->xr = bnow->first-k+1;
bnow->conlist->frac = (bnow->pl[k]-bnow->conlist->pl)/(bnow->pl[k]-bnow->pl[k-1]);
}
}
return;
}
/**********************************************************
Function to assign current
**************************************************************/
void Assign_Current(branch *bnow, double *x, double fac )
{
if (bnow->child != NULL ) Assign_Current(bnow->child, x, fac );
if (bnow->peer != NULL ) Assign_Current(bnow->peer, x, fac );
if ( bnow->conlist != NULL ) {
x[bnow->conlist->xl] -= fac*(bnow->conlist->frac)*(bnow->conlist->amp);
x[bnow->conlist->xr] -= fac*(1.0-(bnow->conlist->frac))*(bnow->conlist->amp);
}
return;
}
/*********************************************************************
Function to count total dendritic length
*********************************************************************/
void count_dendritic_length ( branch *b, double *length )
{
(*length) += b->p_len;
if ( b->child != NULL ) count_dendritic_length ( b->child, length );
if ( b->peer != NULL ) count_dendritic_length ( b->peer, length );
return;
}
/********************************************************************
Function to re-compute physical lengths
********************************************************************/
void phys_lengths( branch *b, double *h )
{
int i,n,k;
double tmp, hnow, *dtmp, *ptmp, *ptr;
n = ((int) ceil((b->p_len)/(*h))) + 1;
hnow = b->p_len/((double) n-1);
dtmp = (double *) malloc( n*sizeof(double));
ptmp = (double *) malloc( n*sizeof(double));
for ( i=0; i<(n-1); i++ ) ptmp[i] = ((double) i)*hnow;
ptmp[n-1] = b->p_len;
dtmp[0] = b->d[0];
for ( i=1 ; i<(n-1) ; i++ ) {
k = 0;
while ( ptmp[i] > b->pl[k] ) k++;
dtmp[i] = b->d[k-1] + (b->d[k] - b->d[k-1])*(ptmp[i] - b->pl[k-1])/(b->pl[k]-b->pl[k-1]);
}
dtmp[n-1] = b->d[b->nobs-1];
ptr = b->d;
b->d = dtmp;
free(ptr);
ptr = b->pl;
b->pl = ptmp;
free(ptr);
b->nobs = n;
if (b->child != NULL) phys_lengths( b->child, h);
if (b->peer != NULL) phys_lengths( b->peer, h);
}
/**********************************************************
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+35.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+60.0)/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+60.0);
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+50.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+60.0);
return 0.125*fac*exp(-tmp);
}
