Skip to content

Commit

Permalink
Merge branch 'stuart-knock-dendrite-differential'
Browse files Browse the repository at this point in the history 
+ This branch solves #43:
    -> Implements a differential form for dendritic dynamics (preexisting RK4)
    -> Leaves the integral form as legacy code.
       -> DendriteIntrgral should be used only for steady state dynamics
          and uniform parameters (alpha and beta)
  • Loading branch information
Paula Sanz-Leon committed Feb 9, 2017
2 parents a7c1af7 + 4e7e2bf commit 903741c
Show file tree
Hide file tree
Showing 28 changed files with 274 additions and 57 deletions.
89 changes: 89 additions & 0 deletions configs/dendrite_integral.conf
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,89 @@
Example configuration file for the Robinson et. al. corticothalamic model

Connection 1 - Excitatory -> Excitatory
Connection 2 - Inhibitory -> Excitatory
Connection 3 - Relay -> Excitatory
Connection 4 - Excitatory -> Inhibitory
Connection 5 - Inhibitory -> Inhibitory
Connection 6 - Relay -> Inhibitory
Connection 7 - Excitatory -> Reticular
Connection 8 - Relay -> Reticular
Connection 9 - Excitatory -> Relay
Connection 10 - Reticular -> Relay
Connection 11 - Stimulus -> Relay

Time: 15 Deltat: 0.0001
Nodes: 144

Connection matrix:
From: 1 2 3 4 5
To 1: 1 2 0 3 0
To 2: 4 5 0 6 0
To 3: 7 0 0 8 0
To 4: 9 0 10 0 11
To 5: 0 0 0 0 0

Population 1: Excitatory
Length: 0.5
Q: 5.248361515
Firing: Sigmoid - Theta: 0.01292 Sigma: 0.0038 Qmax: 340
Dendrite 1: Integral - alpha: 83.33333333 beta: 769.2307692
Dendrite 2: Integral - alpha: 83.33333333 beta: 769.2307692
Dendrite 3: Integral - alpha: 83.33333333 beta: 769.2307692

Population 2: Inhibitory
Length: 0.5
Q: 5.248361515
Firing: Sigmoid - Theta: 0.01292 Sigma: 0.0038 Qmax: 340
Dendrite 4: Integral - alpha: 83.33333333 beta: 769.2307692
Dendrite 5: Integral - alpha: 83.33333333 beta: 769.2307692
Dendrite 6: Integral - alpha: 83.33333333 beta: 769.2307692

Population 3: Reticular
Length: 0.5
Q: 15.39601978
Firing: Sigmoid - Theta: 0.01292 Sigma: 0.0038 Qmax: 340
Dendrite 7: Integral - alpha: 83.33333333 beta: 769.2307692
Dendrite 8: Integral - alpha: 83.33333333 beta: 769.2307692

Population 4: Relay
Length: 0.5
Q: 8.789733431
Firing: Sigmoid - Theta: 0.01292 Sigma: 0.0038 Qmax: 340
Dendrite 9: Integral - alpha: 83.33333333 beta: 769.2307692
Dendrite 10: Integral - alpha: 83.33333333 beta: 769.2307692
Dendrite 11: Integral - alpha: 83.33333333 beta: 769.2307692

Population 5: Stimulation
Length: 0.5
Stimulus: White - Onset: 0 Mean: 1 Psd: 1e-05

Propagator 1: Wave - Tau: 0 Range: 0.086 gamma: 116
Propagator 2: Map - Tau: 0
Propagator 3: Map - Tau: 0.0425
Propagator 4: Wave - Tau: 0 Range: 0.086 gamma: 116
Propagator 5: Map - Tau: 0
Propagator 6: Map - Tau: 0.0425
Propagator 7: Wave - Tau: 0.0425 Range: 0.086 gamma: 116
Propagator 8: Map - Tau: 0
Propagator 9: Wave - Tau: 0.0425 Range: 0.086 gamma: 116
Propagator 10: Map - Tau: 0
Propagator 11: Map - Tau: 0

Coupling 1: Map - nu: 0.001525377176
Coupling 2: Map - nu: -0.003022754434
Coupling 3: Map - nu: 0.0005674779589
Coupling 4: Map - nu: 0.001525377176
Coupling 5: Map - nu: -0.003022754434
Coupling 6: Map - nu: 0.0005674779589
Coupling 7: Map - nu: 0.0001695899041
Coupling 8: Map - nu: 5.070036187e-05
Coupling 9: Map - nu: 0.003447358203
Coupling 10: Map - nu: -0.001465128967
Coupling 11: Map - nu: 0.003593330094

Output: Node: All Start: 5 Interval: 0.5e-2
Population: 1
Dendrite:
Propagator: 1
Coupling:
9 changes: 5 additions & 4 deletions src/de.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,14 +88,15 @@ class RK4 : public Integrator {
void operator=(RK4&);
protected:
double h6; ///< == deltat/6
double deltat5; ///< == deltat*0.5

vector<double> k1;
vector<double> k2;
vector<double> k3;
vector<double> k4;
vector<double> temp;
public:
explicit RK4( DE& de ) : Integrator(de), h6(de.deltat/6.),
explicit RK4( DE& de ) : Integrator(de), h6(de.deltat/6.), deltat5(de.deltat*0.5),
k1(de.n), k2(de.n), k3(de.n), k4(de.n), temp(de.n) {}
~RK4(void) override = default;

Expand All @@ -106,19 +107,19 @@ class RK4 : public Integrator {
}
de.rhs(temp,k1);
for( vvd_size_type i=0; i<de.n; i++ ) {
temp[i] = de.variables[i][j] +de.deltat*.5*k1[i];
temp[i] = de.variables[i][j] +deltat5*k1[i];
}
de.rhs(temp,k2);
for( vvd_size_type i=0; i<de.n; i++ ) {
temp[i] = de.variables[i][j] +de.deltat*.5*k2[i];
temp[i] = de.variables[i][j] +deltat5*k2[i];
}
de.rhs(temp,k3);
for( vvd_size_type i=0; i<de.n; i++ ) {
temp[i] = de.variables[i][j] +de.deltat*k3[i];
}
de.rhs(temp,k4);
for( vvd_size_type i=0; i<de.n; i++ ) {
de.variables[i][j] += h6*( k1[i] +2*k2[i] +2*k3[i] +k4[i] );
de.variables[i][j] += h6*( k1[i] +2.0*k2[i] +2.0*k3[i] +k4[i] );
}
}
}
Expand Down
82 changes: 52 additions & 30 deletions src/dendrite.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,26 @@
#include "dendrite.h"
using std::endl;

/**
@brief Computes the derivatives of the dendritic response function.
The dendritic response function is given by:
\f{eqnarray*}{
\frac{dv}{dt}&=&W \\
\frac{dW}{dt}&=&\left(\nu\phi - V - \left(\frac{1}{\alpha} + \frac{1}{\beta}\right) W\right) \alpha \beta \\
\frac{d\nu\phi}{dt}&=&0
\f}
*/
void Dendrite::DendriteDE::rhs( const vector<double>& y, vector<double>& dydt ) {
// y = {V,W==dv/dt,nuphi}
// dydt = {dv/dt==W, dW/dt==d^2V/dt^2,dnuphi/dt} d(nuphi)/dt from precouple
dydt[0] = y[1];

dydt[1] = (y[2] - y[0] - factorab * y[1]) * alphaxbeta;

dydt[2] = 0.0;
}

void Dendrite::init( Configf& configf ) {
if( !configf.next( label("Dendrite ",index+1) )) {
std::cerr<<"Dendrite "<<index+1<<" not found."<<endl;
Expand All @@ -23,51 +43,53 @@ void Dendrite::init( Configf& configf ) {
} else {
v.resize(nodes,atof(buffer.c_str()));
}
oldnp = v;

configf.param("alpha",alpha);
configf.param("beta",beta);

aminusb = alpha - beta;
expa = exp(-alpha*deltat);
expb = exp(-beta*deltat);
// Initialize constant factors to speed up computation.
factorab = 1./alpha + 1./beta;
alphaxbeta = alpha * beta;

de->init(v[0]); // call Dendrite::DendriteDE::init
de->factorab = factorab;
de->alphaxbeta = alphaxbeta;
}

void Dendrite::DendriteDE::init(const double vinit) {
variables[0].clear();
variables[1].clear();
variables[2].clear();
variables[0].resize(nodes, vinit);
variables[1].resize(nodes, 0.0);
variables[2].resize(nodes, 0.0);
}

Dendrite::Dendrite( size_type nodes, double deltat, size_type index,
const Propagator& prepropag, const Coupling& precouple )
: NF(nodes,deltat,index), v(nodes), dvdt(nodes,0), oldnp(nodes),
prepropag(prepropag), precouple(precouple) {
: NF(nodes,deltat,index), v(nodes), prepropag(prepropag), precouple(precouple) {
de = new DendriteDE(nodes, deltat);
rk4 = new RK4(*de);
}

Dendrite::~Dendrite() = default;
Dendrite::~Dendrite() {
delete de;
delete rk4;
}

void Dendrite::step() {
// assume that alpha, beta are constant and nu*phi is linear for the time step
if(alpha!=beta) {
for(size_type i=0; i<nodes; i++) {
dpdt = ( precouple[i] -oldnp[i] )/deltat;
adjustednp = oldnp[i] -factorab*dpdt -v[i];
C1 = ( adjustednp*beta -dvdt[i] +dpdt )/aminusb;
C1expa = C1*expa;
C2expb = expb*(-C1-adjustednp);
v[i] = C1expa+C2expb+precouple[i] -factorab*dpdt;
dvdt[i] = C1expa*(-alpha) +C2expb*(-beta)+dpdt;
oldnp[i]=precouple[i]; //Save current pulse density for next step
}
} else { // alpha==beta
for(size_type i=0; i<nodes; i++) {
dpdt = ( precouple[i] -oldnp[i] )/deltat;
adjustednp = oldnp[i] -factorab*dpdt -v[i];
C1 = dvdt[i] -alpha*adjustednp -dpdt;
C1dtplusC2 = C1*deltat -adjustednp;
v[i] = C1dtplusC2*expa +precouple[i] -factorab*dpdt;
dvdt[i] = (C1-alpha*C1dtplusC2)*expa +dpdt;
oldnp[i]=precouple[i]; //Save current pulse density for next step
}
//Copy the \nu\phi into their corresponding variable in the DE.
for( size_type i=0; i<nodes; i++ ) {
(*de)[2][i] = precouple[i]; // \nu\phi
}
//Integrate the dendritic response one step forward in time.
rk4->step();
//Copy the voltages from the updated DE to the local variable v.
for( size_type i=0; i<nodes; i++ ) {
v[i] = (*de)[0][i]; // Voltage
}
}


void Dendrite::output( Output& output ) const {
output("Dendrite",index+1,"V",v);
}
40 changes: 18 additions & 22 deletions src/dendrite.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,29 +17,25 @@ class Dendrite : public NF {
Dendrite(void); // default constructor
Dendrite(Dendrite& ); // no copy constructor

// variables that are intialized once to speed up computation
double aminusb; ///< == alpha - beta
double expa; ///< == exp(-alpha*deltat)
double expb; ///< == exp(-beta*deltat)
double factorab; ///< == 1./alpha + 1./beta;

// variables that are used every timestep
double adjustednp;
double deltaPdeltat;
double C1;
double dpdt;
double C1expa;
double C2expb;
double C1dtplusC2;
protected:

double alpha; // needed here for DendriteRamp
double beta;
// variables that are initialized once to speed up computation
double factorab; ///< == 1./alpha + 1./beta;
double alphaxbeta; ///< == alpha * beta;

vector<double> v;
vector<double> dvdt;
//vector<double> np;
vector<double> oldnp;
protected:
struct DendriteDE : public DE {
double factorab, alphaxbeta;
virtual void init( const double vinit);
DendriteDE( size_type nodes, double deltat) : DE(nodes, deltat, 3) {}
~DendriteDE(void) override = default;
void rhs( const vector<double>& y, vector<double>& dydt ) override;
};
DendriteDE* de;
RK4* rk4;

double alpha; ///< Mean decay rate of the soma response to a delta-function synaptic input (needed here for DendriteRamp).
double beta; ///< Mean rise rate of the soma response to a delta-function synaptic input.

vector<double> v; ///< Membrane potential.

void init( Configf& configf ) override;
//virtual void restart( Restartf& restartf );
Expand Down
56 changes: 56 additions & 0 deletions src/dendrite_integral.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
/** @file dendrite_integral.cpp
@brief DendriteIntegral handles the dendritic response of the postsynaptic
population using the integral form for the step() member-function.
A more detailed multiline description...
@author Peter Drysdale, Felix Fung,
*/

#include "dendrite_integral.h"


void DendriteIntegral::init( Configf& configf ) {
Dendrite::init(configf);

oldnp = v;

// Initialize constant factors to speed up computation.
aminusb = alpha - beta;
expa = exp(-alpha*deltat);
expb = exp(-beta*deltat);
factorab = 1./alpha + 1./beta;
}

DendriteIntegral::DendriteIntegral( size_type nodes, double deltat, size_type index,
const Propagator& prepropag, const Coupling& precouple )
: Dendrite(nodes, deltat, index, prepropag, precouple), dvdt(nodes,0), oldnp(nodes) {
}

DendriteIntegral::~DendriteIntegral() = default;

void DendriteIntegral::step() {
// assume that alpha, beta are constant and nu*phi is linear for the time step
if(alpha!=beta) {
for(size_type i=0; i<nodes; i++) {
dpdt = ( precouple[i] -oldnp[i] )/deltat;
adjustednp = oldnp[i] -factorab*dpdt -v[i];
C1 = ( adjustednp*beta -dvdt[i] +dpdt )/aminusb;
C1expa = C1*expa;
C2expb = expb*(-C1-adjustednp);
v[i] = C1expa+C2expb+precouple[i] -factorab*dpdt;
dvdt[i] = C1expa*(-alpha) +C2expb*(-beta)+dpdt;
oldnp[i]=precouple[i]; //Save current pulse density for next step
}
} else { // alpha==beta
for(size_type i=0; i<nodes; i++) {
dpdt = ( precouple[i] -oldnp[i] )/deltat;
adjustednp = oldnp[i] -factorab*dpdt -v[i];
C1 = dvdt[i] -alpha*adjustednp -dpdt;
C1dtplusC2 = C1*deltat -adjustednp;
v[i] = C1dtplusC2*expa +precouple[i] -factorab*dpdt;
dvdt[i] = (C1-alpha*C1dtplusC2)*expa +dpdt;
oldnp[i]=precouple[i]; //Save current pulse density for next step
}
}
}
50 changes: 50 additions & 0 deletions src/dendrite_integral.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@
/** @file dendrite_integral.h
@brief Defines the DendriteIntegral class, which handles the response of the
postsynaptic population using the integral form for the step()
member-function.
A more detailed multiline description...
@author Peter Drysdale, Felix Fung,
*/

#ifndef NEUROFIELD_SRC_DENDRITE_INTEGRAL_H
#define NEUROFIELD_SRC_DENDRITE_INTEGRAL_H

#include "dendrite.h"

class DendriteIntegral : public Dendrite {
DendriteIntegral(void); // default constructor
DendriteIntegral(DendriteIntegral& ); // no copy constructor

// variables that are intialized once to speed up computation
double aminusb; ///< == alpha - beta
double expa; ///< == exp(-alpha*deltat)
double expb; ///< == exp(-beta*deltat)
double factorab; ///< == 1./alpha + 1./beta;

// variables that are used every timestep
double adjustednp;
double deltaPdeltat; //NOTE: doesn't seem to be used anywhere (???maybe replaced by current dpdt???).
double C1;
double dpdt; // Temporal derivative of pulse density.
double C1expa;
double C2expb;
double C1dtplusC2;

protected:
vector<double> dvdt; ///< Temporal derivative of the membrane potential.
//vector<double> np; ///< Pulse density (nu*phi).
vector<double> oldnp; ///< Pulse density (nu*phi) of the previous time-step.

void init( Configf& configf ) override;

public:

DendriteIntegral( size_type nodes, double deltat, size_type index,
const Propagator& prepropag, const Coupling& precouple );
~DendriteIntegral(void) override;
void step(void) override;
};

#endif //NEUROFIELD_SRC_DENDRITE_INTEGRAL_H
Loading

0 comments on commit 903741c

Please sign in to comment.