SHDynamics.m

Go to the documentation of this file.

namespace models{
namespace motorunit{


/* (Autoinserted by mtoc++)
 * This source code has been filtered by the mtoc++ executable,
 * which generates code that can be processed by the doxygen documentation tool.
 *
 * On the other hand, it can neither be interpreted by MATLAB, nor can it be compiled with a C++ compiler.
 * Except for the comments, the function bodies of your M-file functions are untouched.
 * Consequently, the FILTER_SOURCE_FILES doxygen switch (default in our Doxyfile.template) will produce
 * attached source files that are highly readable by humans.
 *
 * Additionally, links in the doxygen generated documentation to the source code of functions and class members refer to
 * the correct locations in the source code browser.
 * However, the line numbers most likely do not correspond to the line numbers in the original MATLAB source files.
 */

class SHDynamics
  :public dscomponents.ACoreFun {
    public:

        double MSLink_MaxFactorSignal = 40;
        double MSLink_MaxFactor = 7;
        double MSLink_MinFactor = .3;
    public:

        SarcoConst_slow;
        SarcoConst_fast;
        SarcoConst_base;
        .-dependent SarcoConst_dynpos;
    private:

        spm;

        
    public: /* ( Transient ) */

        hadpeak = false;

        motoconst;

        sarcoconst;

        signalweight;

        
    public: /* ( Transient ) */

        SHDynamics(sys) {
            this = this@dscomponents.ACoreFun(sys);

            this.initSarcoConst;
            this.initJSparsityPattern;
            this.fDim= sys.dsa + sys.dm;
            this.xDim= this.fDim;
            this.spm= sys.Model.SinglePeakMode;
        }


        function  prepareSimulation(colvec<double> mu) {
            prepareSimulation@dscomponents.ACoreFun(this, mu);

            /*  Precompute motoneuron/sarcomere constants */
            this.motoconst= this.getMotoConst(mu(1));
            this.sarcoconst= this.getSarcoConst(mu(1));

            this.hadpeak= false;
        }


        function dy = evaluate(matrix<double>, y,rowvec<double> t) {

            dm = this.System.dm;

            /* % Inner dynamics (with partial nonlinear linking) */
            dy = zeros(size(y));

            /* % Neuro dynamics (last 2 arguments are 2 phases: phase_prim, phase_sec) */
            c = this.motoconst;
            /*  dendrites */
            dy(1) = (-c(1)*(y(1)-c(11))-c(5)*(y(1)-y(2)))/c(7);
            /*  soma */
            dy(2) = (-c(6)*(y(2)-c(11))-c(5)*(y(2)-y(1))...
                -c(4)*y(3).^3*y(4)*(y(2)-c(9))...
                -c(2)*y(5).^4*(y(2)-c(10))...
                -c(3)*y(6).^2*(y(2)-c(10)))/c(8);
            /*  the four gating variables */
            dy(3) = 0.32*(13-y(2))/(exp((13-y(2))/5)-1)*(1-y(3))...
                -0.28*(y(2)-40)/(exp((y(2)-40)/5)-1)*(y(3));
            dy(4) = 0.128*(exp((17-y(2))/18))*(1-y(4))-4/(exp((40-y(2))/5)+1)*(y(4));
            dy(5) = 0.032*(15-y(2))/(exp((15-y(2))/5)-1)*(1-y(5))...
                -0.5*(exp((10-y(2))/40))*(y(5));
            dy(6) = 3.5/(exp((55-y(2))/4)+1)*(1-y(6))-0.025*(y(6));

            /* % Sarcomere dynamics */
            dy(dm+1:end) = this.dydt_sarcomere(y(dm+1:end),t);

            /* % Link of motoneuron to sarcomer cell */
            if this.spm && this.hadpeak
                signal = 0;
            else
                fac = this.MSLink_MaxFactor;
                if y(2) < this.MSLink_MaxFactorSignal
                    fac = this.getLinkFactor(y(2));
                end
                /*  "Link" at the only sarcomere */
                signal = fac*y(2)/this.sarcoconst(1);
            end

            dy(dm+1,:) = dy(dm+1,:) + signal;
        }
        function dy = evaluateCoreFun(matrix<double> y,double t,colvec<double> mu) {
            error(" evaluate is overridden directly. ");
        }


        function J = getStateJacobian(matrix<double>, y,rowvec<double> t) {

            s = this.System;
            dm = s.dm;
            dsa = s.dsa;
            bs = dm + dsa;
            J = spalloc(bs,bs,45+303);  /*  45: 40 entries in motoneuron and spindle part + reserve, 303: entries per single sarcomere */


            /* % Neuron jacobian */
            c = this.motoconst;
            JN = zeros(6,6);
            JN(1,1) = -(c(1) + c(5))/c(7);
            JN(2,1) = c(5)/c(8);
            JN(1,2) = c(5)/c(7);
            JN(2,2) = -(c(4)*y(4)*y(3)^3 + c(2)*y(5)^4 + c(3)*y(6)^2 + c(5) + c(6))/c(8);
            JN(3,2) = (8*(y(3) - 1))/(25*(exp(13/5 - y(2)/5) - 1)) - (7*y(3))/(25*(exp(y(2)/5 - 8) - 1)) + (y(3)*exp(y(2)/5 - 8)*((7*y(2))/25 - 56/5))/(5*(exp(y(2)/5 - 8) - 1)^2) + (exp(13/5 - y(2)/5)*((8*y(2))/25 - 104/25)*(y(3) - 1))/(5*(exp(13/5 - y(2)/5) - 1)^2);
            JN(4,2) = (8*exp(17/18 - y(2)/18)*(y(4) - 1))/1125 - (4*y(4)*exp(8 - y(2)/5))/(5*(exp(8 - y(2)/5) + 1)^2);
            JN(5,2) = (y(5)*exp(1/4 - y(2)/40))/80 + (4*(y(5) - 1))/(125*(exp(3 - y(2)/5) - 1)) + (exp(3 - y(2)/5)*((4*y(2))/125 - 12/25)*(y(5) - 1))/(5*(exp(3 - y(2)/5) - 1)^2);
            JN(6,2) = -(7*exp(55/4 - y(2)/4)*(y(6) - 1))/(8*(exp(55/4 - y(2)/4) + 1)^2);
            JN(2,3) = (3*c(4)*y(3)^2*y(4)*(c(9) - y(2)))/c(8);
            JN(3,3) =  ((8*y(2))/25 - 104/25)/(exp(13/5 - y(2)/5) - 1) - ((7*y(2))/25 - 56/5)/(exp(y(2)/5 - 8) - 1);
            JN(2,4) = (c(4)*y(3)^3*(c(9) - y(2)))/c(8);
            JN(4,4) =  - (16*exp(17/18 - y(2)/18))/125 - 4/(exp(8 - y(2)/5) + 1);
            JN(2,5) = (4*c(2)*y(5)^3*(c(10) - y(2)))/c(8);
            JN(5,5) = ((4*y(2))/125 - 12/25)/(exp(3 - y(2)/5) - 1) - exp(1/4 - y(2)/40)/2;
            JN(2,6) = (2*c(3)*y(6)*(c(10) - y(2)))/c(8);
            JN(6,6) = - 7/(2*(exp(55/4 - y(2)/4) + 1)) - 1/40;
            J(1:6,1:6) = JN;

            /* % Sarcomere jacobian */
            J(dm+1:bs,dm+1:bs) = this.Jac_Sarco(y(dm+1:end), t);
            fac = this.MSLink_MaxFactor;
            if (y(2) < this.MSLink_MaxFactorSignal)
                fac = this.getLinkFactor(y(2));
            end
            /*  link position is dm+1 */
            J(dm+1,2) = fac/this.sarcoconst(1);
        }
        function copy = clone() {
            copy = models.motorunit.Dynamics(this.System);

            /*  Call superclass clone (for deep copy) */
            copy = clone@dscomponents.ACompEvalCoreFun(this, copy);

            /*  Copy local properties
             * No local properties are to be copied here, as so far
             * everything is done in the constructor. */
        }
        function  initJSparsityPattern() {
            s = this.System;
            dm = s.dm; dsa = s.dsa;
            size = dm + dsa;
            J = sparse(size,size);

            /*  Neuro */
            i = [1,1,2,2,2,2,2,2,3,3,4,4,5,5,6,6];
            j = [1,2,1,2,3,4,5,6,2,3,2,4,2,5,2,6];
            J(1:dm,1:dm) = sparse(i,j,true,6,6);

            /*  Sarco */
            mc = metaclass(this);
            s = load(fullfile(fileparts(which(mc.Name))," JSparsityPattern "));
            J(dm+1:size,dm+1:size) = s.JP;

            linkpos = dm+1;
            J(linkpos,2) = true;
            this.JSparsityPattern= logical(J);
        }
        function f = getLinkFactor(rowvec<double> y) {
            f = this.MSLink_MinFactor + exp(-(y-this.MSLink_MaxFactorSignal).^2/150)...
                *(this.MSLink_MaxFactor-this.MSLink_MinFactor);
        }
    private: /* ( Transient ) */

        function sc = getSarcoConst(mu_fibretype) {
            sc = this.SarcoConst_base;
            sc(this.SarcoConst_dynpos) = (1-mu_fibretype(1))*this.SarcoConst_slow + mu_fibretype(1)*this.SarcoConst_fast;
        }
        function c = getMotoConst(mu_fibretype) {

            /*  Membrane capacitance (NOT to be confused with Cm of the
             * sarcomere model!) */
            Cm=1;

            Ri=70/1000;
            c = zeros(11,length(mu_fibretype));
            /*  cf. Cisi and Kohn 2008, Table 2, page 7 */
            Rmd = 14.4+6.05-coolExp(6.05,14.4,mu_fibretype);
            Rms=1.15+0.65-coolExp(0.65,1.15,mu_fibretype);

            ld=coolExp(0.55,1.06,mu_fibretype);
            ls=coolExp(77.5e-6*100,113e-6*100,mu_fibretype);

            rd=coolExp(41.5e-6*100,92.5e-6*100,mu_fibretype)/2;
            rs=coolExp(77.5e-6*100,113e-6*100,mu_fibretype)/2;

            c(1,:) = 2*pi*rd.*ld./Rmd;   /*  para.Gld */

            c(2,:) = 4*2*pi*rs.*ls;      /*  para.Gkf */

            c(3,:) = 16*2*pi*rs.*ls;     /*  para.Gks */

            c(4,:) = 30*2*pi*rs.*ls;     /*  para.Gna */

            c(5,:) = 2./(Ri*ld./(pi*rd.^2)+Ri*ls./(pi*rs.^2));     /*  para.Gc */

            c(6,:) = 2*pi*rs.*ls./Rms;   /*  para.Gls */

            c(7,:) = 2*pi*rd.*ld*Cm;     /*  para.Cd */

            c(8,:) = 2*pi*rs.*ls*Cm;     /*  para.Cs */

            s = ones(size(mu_fibretype));
            c(9,:) = 120*s;     /*  para.Vna */

            c(10,:) = -10*s;     /*  para.Vk */

            c(11,:) = 0*s;     /*  para.Vl */


            function v = coolExp(a,b,mu)
                v = exp(log(100)*mu)*(b-a)/100 + a;
            end
        }
        function dy = dydt_sarcomere(matrix<double> y,rowvec<double> t);

        function  initSarcoConst();

        function J = Jac_Sarco(matrix<double> y,rowvec<double> t);

};
}
}