diff --git a/configs/dendrite_integral.conf b/configs/dendrite_integral.conf new file mode 100644 index 0000000..822d05a --- /dev/null +++ b/configs/dendrite_integral.conf @@ -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: diff --git a/src/de.h b/src/de.h index 3c8f1aa..a963790 100644 --- a/src/de.h +++ b/src/de.h @@ -88,6 +88,7 @@ class RK4 : public Integrator { void operator=(RK4&); protected: double h6; ///< == deltat/6 + double deltat5; ///< == deltat*0.5 vector k1; vector k2; @@ -95,7 +96,7 @@ class RK4 : public Integrator { vector k4; vector 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; @@ -106,11 +107,11 @@ class RK4 : public Integrator { } de.rhs(temp,k1); for( vvd_size_type i=0; i& y, vector& 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 "<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; istep(); + //Copy the voltages from the updated DE to the local variable v. + for( size_type i=0; i v; - vector dvdt; - //vector np; - vector 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& y, vector& 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 v; ///< Membrane potential. void init( Configf& configf ) override; //virtual void restart( Restartf& restartf ); diff --git a/src/dendrite_integral.cpp b/src/dendrite_integral.cpp new file mode 100644 index 0000000..1da2d5f --- /dev/null +++ b/src/dendrite_integral.cpp @@ -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 dvdt; ///< Temporal derivative of the membrane potential. + //vector np; ///< Pulse density (nu*phi). + vector 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 diff --git a/src/dendriteramp.h b/src/dendriteramp.h index 7458ac5..2d4119a 100644 --- a/src/dendriteramp.h +++ b/src/dendriteramp.h @@ -40,7 +40,7 @@ class DendriteRamp : public Dendrite { void step(void) override; inline const vector& V(void) const { - return (*de)[1]; + return (*de)[1]; //NOTE: Pretty sure this is an error: index=1 points to dv/dt, to get voltage, which is what I think is wanted here, it should be index=0 } void output( Output& output ) const override; }; diff --git a/src/qresponse.cpp b/src/qresponse.cpp index 0c9ba82..9fc9267 100644 --- a/src/qresponse.cpp +++ b/src/qresponse.cpp @@ -15,6 +15,7 @@ #include "fmath.h" #include "dendriteramp.h" +#include "dendrite_integral.h" using std::endl; void QResponse::init( Configf& configf ) { @@ -106,6 +107,8 @@ void QResponse::add2Dendrite( size_type index, // PUT YOUR DENDRITE HERE if( temp == "Ramp" ) { dendrites.add( new DendriteRamp(nodes,deltat,index,prepropag,precouple) ); + } else if ( temp == "Integral" ) { + dendrites.add( new DendriteIntegral(nodes,deltat,index,prepropag,precouple) ); // END PUT YOUR DENDRITE HERE } else { dendrites.add( new Dendrite(nodes,deltat,index,prepropag,precouple) ); diff --git a/test/data/configs/Plasticity/bcm.output.gz b/test/data/configs/Plasticity/bcm.output.gz index b86f05e..037c25c 100644 Binary files a/test/data/configs/Plasticity/bcm.output.gz and b/test/data/configs/Plasticity/bcm.output.gz differ diff --git a/test/data/configs/Plasticity/bihemispheric.output.gz b/test/data/configs/Plasticity/bihemispheric.output.gz index 6fc0273..0b38037 100644 Binary files a/test/data/configs/Plasticity/bihemispheric.output.gz and b/test/data/configs/Plasticity/bihemispheric.output.gz differ diff --git a/test/data/configs/Plasticity/eirs-bcm.output.gz b/test/data/configs/Plasticity/eirs-bcm.output.gz index 41300f3..56edd0f 100644 Binary files a/test/data/configs/Plasticity/eirs-bcm.output.gz and b/test/data/configs/Plasticity/eirs-bcm.output.gz differ diff --git a/test/data/configs/Plasticity/marcus.output.gz b/test/data/configs/Plasticity/marcus.output.gz index 6f4b1b1..0908e43 100644 Binary files a/test/data/configs/Plasticity/marcus.output.gz and b/test/data/configs/Plasticity/marcus.output.gz differ diff --git a/test/data/configs/Plasticity/network.output.gz b/test/data/configs/Plasticity/network.output.gz index eb44f5e..35db7e7 100644 Binary files a/test/data/configs/Plasticity/network.output.gz and b/test/data/configs/Plasticity/network.output.gz differ diff --git a/test/data/configs/Plasticity/plasticity.output.gz b/test/data/configs/Plasticity/plasticity.output.gz index 7e56c5e..998aa35 100644 Binary files a/test/data/configs/Plasticity/plasticity.output.gz and b/test/data/configs/Plasticity/plasticity.output.gz differ diff --git a/test/data/configs/Plasticity/rtms.output.gz b/test/data/configs/Plasticity/rtms.output.gz index a1c19ae..dcd0f2e 100644 Binary files a/test/data/configs/Plasticity/rtms.output.gz and b/test/data/configs/Plasticity/rtms.output.gz differ diff --git a/test/data/configs/Plasticity/tbs.output.gz b/test/data/configs/Plasticity/tbs.output.gz index 26773c1..b50df02 100644 Binary files a/test/data/configs/Plasticity/tbs.output.gz and b/test/data/configs/Plasticity/tbs.output.gz differ diff --git a/test/data/configs/TMS/teps-bi.output.gz b/test/data/configs/TMS/teps-bi.output.gz index 51d257b..abd39b7 100644 Binary files a/test/data/configs/TMS/teps-bi.output.gz and b/test/data/configs/TMS/teps-bi.output.gz differ diff --git a/test/data/configs/TMS/teps-gabab.output.gz b/test/data/configs/TMS/teps-gabab.output.gz index c2c02db..2512bef 100644 Binary files a/test/data/configs/TMS/teps-gabab.output.gz and b/test/data/configs/TMS/teps-gabab.output.gz differ diff --git a/test/data/configs/TMS/teps-mono.output.gz b/test/data/configs/TMS/teps-mono.output.gz index 6eda0da..07757e1 100644 Binary files a/test/data/configs/TMS/teps-mono.output.gz and b/test/data/configs/TMS/teps-mono.output.gz differ diff --git a/test/data/configs/aerp.output.gz b/test/data/configs/aerp.output.gz index ae9197f..6d61aab 100644 Binary files a/test/data/configs/aerp.output.gz and b/test/data/configs/aerp.output.gz differ diff --git a/test/data/configs/burst.output.gz b/test/data/configs/burst.output.gz index 927fd89..b591e16 100644 Binary files a/test/data/configs/burst.output.gz and b/test/data/configs/burst.output.gz differ diff --git a/test/data/configs/dendrite_integral.output.gz b/test/data/configs/dendrite_integral.output.gz new file mode 100644 index 0000000..4f49f91 Binary files /dev/null and b/test/data/configs/dendrite_integral.output.gz differ diff --git a/test/data/configs/diff_arctan.output.gz b/test/data/configs/diff_arctan.output.gz index 6344883..8e0220a 100644 Binary files a/test/data/configs/diff_arctan.output.gz and b/test/data/configs/diff_arctan.output.gz differ diff --git a/test/data/configs/emirs.output.gz b/test/data/configs/emirs.output.gz index 7feecdb..f3f0c4d 100644 Binary files a/test/data/configs/emirs.output.gz and b/test/data/configs/emirs.output.gz differ diff --git a/test/data/configs/example.output.gz b/test/data/configs/example.output.gz index ef94c65..d7085dc 100644 Binary files a/test/data/configs/example.output.gz and b/test/data/configs/example.output.gz differ diff --git a/test/data/configs/onepop.output.gz b/test/data/configs/onepop.output.gz index 997060e..28a2be8 100644 Binary files a/test/data/configs/onepop.output.gz and b/test/data/configs/onepop.output.gz differ diff --git a/test/data/configs/ramp.output.gz b/test/data/configs/ramp.output.gz index 485f17c..584b480 100644 Binary files a/test/data/configs/ramp.output.gz and b/test/data/configs/ramp.output.gz differ diff --git a/test/data/configs/wave.output.gz b/test/data/configs/wave.output.gz index 596f60d..d282dc0 100644 Binary files a/test/data/configs/wave.output.gz and b/test/data/configs/wave.output.gz differ