membraneModels.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2013-2016 Zuse Institute Berlin                            */
/*  Copyright (C) 2013-2013 U Milan                                          */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef MEMBRANEMODELS_HH
#define MEMBRANEMODELS_HH
 
#include <dune/common/config.h>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
 
 
namespace Kaskade {
  
  template <class Derived, int n>
  struct MembraneModelBase
  {
    static int const nGating = n;
    
    typedef Dune::FieldVector<double,nGating>         Gating;
    
    typedef Dune::FieldMatrix<double,nGating,nGating> GatingJacobian;
    
    MembraneModelBase() {}
    
    MembraneModelBase(std::string const& name): nm(name) {}
 
    //MembraneModelBase(std::string const& name, double uFix, Gating const& vFix): nm(name), fix(uFix,vFix) {}
    
    std::string const& name() const { return nm; }
    
    std::pair<double,Gating> const& restState() const { return fix; }
    
    double current_du(double u, Gating const& v, double h=1e-5) const
    {
      double fuv = static_cast<Derived const&>(*this).current(u-h,v);
      double fuiv = static_cast<Derived const&>(*this).current(u+h,v);
      
      return (fuiv-fuv)/(2*h);
    }
    
    Gating current_dv(double u, Gating const& v) const
    {
      double f = static_cast<Derived const&>(*this).current(u,v);
      Gating fv_uv,dir;
      for (int i=0; i<nGating; ++i)
      {
        dir = 0;
        dir[i] = 1e-5;
        fv_uv[i] = (static_cast<Derived const&>(*this).current(u,v+dir)-f)/1e-5;
      }
      return fv_uv;
    }
        
    Gating gatingRhs_du(double u, Gating const& v) const
    {
      Gating guA = static_cast<Derived const&>(*this).gatingRhs(u-1e-5,v);
      Gating guB = static_cast<Derived const&>(*this).gatingRhs(u+1e-5,v);
      return (guB-guA)/2e-5;
    }
    
    GatingJacobian gatingRhs_dv(double u, Gating const& v) const
    {
      GatingJacobian gv_uv;
      Gating col, dir, guv = static_cast<Derived const&>(*this).gatingRhs(u,v);
      for (int i=0; i<nGating; ++i){
        dir = 0;
        dir[i]=1e-5;
        col = (static_cast<Derived const&>(*this).gatingRhs(u,v+dir)-guv)/1e-5;
        for (int j=0; j<nGating; ++j) gv_uv[j][i] = col[j];
      }
      return gv_uv;
    }
 
  protected:
    void setRestState(double uFix, Gating const& vFix)
    {
      using State = Dune::FieldVector<double,nGating+1>;
      using DState = Dune::FieldMatrix<double,nGating+1,nGating+1>;
 
      // Get the initial value for Newton's method
      State x, dx; 
      x[0] = uFix; for (int i=0; i<nGating; ++i) x[i+1] = vFix[i];
      std::cout << "x_init = " << x << std::endl;
      
      State F; F = 0;
      DState DF; DF = 0;
      
      // perform a few ordinary (or simplified?) Newton iterations for the system [current,gatingRhs] = 0.    
    
      double uFix0 = uFix;
      Gating vFix0(vFix);
      
      double normF = 0, dxNorm = 0;
      int steps = 0;
      
      do
      {
        steps++;
        
        double c = static_cast<Derived const&>(*this).current(uFix0,vFix0);
      
        Gating g = static_cast<Derived const&>(*this).gatingRhs(uFix0,vFix0);
       
        F[0] = c; for (int i=0; i<nGating; ++i) F[i+1] = g[i];
        normF = F.infinity_norm();
 
        double current_u = current_du(uFix0, vFix0);
        Gating current_v = static_cast<Derived const&>(*this).current_dv(uFix0, vFix0);
        Gating         gatingRHS_u = static_cast<Derived const&>(*this).gatingRhs_du(uFix0, vFix0);
        GatingJacobian gatingRHS_v = gatingRhs_dv(uFix0, vFix0);
        
       
        DF[0][0] = current_u;
        for (int i=0; i<nGating; ++i) DF[0][i+1] = -current_v[i];
        for (int i=0; i<nGating; ++i) DF[i+1][0] =  gatingRHS_u[i];   
        for (int i=0; i<nGating; ++i) 
          for (int j=0; j<nGating; ++j) 
            DF[i+1][j+1] = gatingRHS_v[i][j];
 
        
        // solve DF*dx = -F, dx = x(j+1) - x(j)
        DF.solve(dx,-F);
        
        double d = 0.1;
        uFix0 += d*dx[0];
        for (int i=0; i<nGating; ++i)
          vFix0[i] += d*dx[i+1];
        
        dxNorm = dx.infinity_norm();
      } while ( (normF > 1e-6) || (dxNorm > 1e-4) );
      
      std::cout << " Newton steps = " << steps << std::endl;
      
      // store the value
      fix.first = uFix0;
      for (int i=0; i<nGating; ++i)
        fix.second[i] = vFix0[i];
                
      std::cout << "fixpoint = " << fix.first << "   " << fix.second << std::endl;    
 
    }
 
    
  private:
    std::string nm;
    std::pair<double,Gating> fix;
  };
  
  //----------------------------------------------------------------------------------------
  
  
  struct AlievPanfilov: public MembraneModelBase<AlievPanfilov,1>
  {
    static int const nGating = 1;
 
    typedef Dune::FieldVector<double,nGating> Gating;
    typedef Dune::FieldMatrix<double,nGating,nGating> GatingJacobian;
 
    AlievPanfilov();
 
    double current(double u, Gating const& v) const;
 
    double current_du(double u, Gating const& v) const;
 
    Gating current_dv(double u, Gating const& v) const;
 
    Gating gatingRhs(double u, Gating const& v, double Istim=0) const;
 
    Gating gatingRhs_du(double u, Gating const& v) const;
 
    GatingJacobian gatingRhs_dv(double u, Gating const& v) const;
 
  private:
    double a, eps0, k, mu1, mu2;
  };
 
  //-----------------------------------------------------------------------------------------------------
 
  struct PhysiologicalAlievPanfilov: public MembraneModelBase<PhysiologicalAlievPanfilov,1>
  {
    static int const nGating = 1;
 
    typedef Dune::FieldVector<double,nGating> Gating;
    typedef Dune::FieldMatrix<double,nGating,nGating> GatingJacobian;
 
    PhysiologicalAlievPanfilov(double Cm = 1e-2/16);
 
    double current(double u, Gating const& v) const;
 
    double current_du(double u, Gating const& v) const;
 
    Gating current_dv(double u, Gating const& v) const;
 
    Gating gatingRhs(double u, Gating const& v, double Istim=0) const;
 
    Gating gatingRhs_du(double u, Gating const& v) const;
 
    GatingJacobian gatingRhs_dv(double u, Gating const& v) const;
 
  private:
    AlievPanfilov ap;
    double Cm, s;
  };
 
  //-----------------------------------------------------------------------------------------------------
 
  class FentonKarma: public MembraneModelBase<FentonKarma,2>
  {
  public:
    static int const nGating = 2;
 
    enum ParameterSet { BR      
                      , MBR     
                      , MLRI    
                      , GP      
                      };
 
    typedef Dune::FieldVector<double,nGating> Gating;
    typedef Dune::FieldMatrix<double,nGating,nGating> GatingJacobian;
 
    FentonKarma(ParameterSet p, double Cm = 1e-2);
 
    double current(double u, Gating const& v) const;
 
    double current_du(double u, Gating const& v) const;
 
    Gating current_dv(double u, Gating const& v) const;
 
    Gating gatingRhs(double u, Gating const& v, double Istim=0) const;
 
    Gating gatingRhs_du(double u, Gating const& v) const;
 
    GatingJacobian gatingRhs_dv(double u, Gating const& v) const;
 
  private:
    // A smoothed Heaviside function. The original Fenton-Karma model is partially formulated
    // in terms of the discontinuous Heaviside function, which is in general a bad idea for numerical
    // computations.
    double h(double x) const;
 
    double u0 = -0.085; // resting voltage
    double ufi = 0.015; // maximum voltage
    double k = 10;      // scaling factor for tanh-based smooth Heaviside approximation as in article
    double K = 3;       // additional scaling factor for approximating non-smooth Heaviside
    double Cm, gfi, taur, tausi, tau0, tauvplus, tauv1minus, tauv2minus, tauwplus, tauwminus, uc, uv, usic, taud;
  };
 
  //-----------------------------------------------------------------------------------------------------
 
  struct TenTusscher16: public MembraneModelBase<TenTusscher16,16>
  {
    TenTusscher16();
    
    double current(double u, Gating const& v) const;
 
    double current_du(double u, Gating const& v, double h=1e-5) const;
    
    Gating current_dv(double u, Gating const& v) const;
    
    Gating gatingRhs(double u, Gating const& v, double Istim=0) const;
    
    Gating gatingRhs_du(double u, Gating const& v) const;
    
    GatingJacobian gatingRhs_dv(double u, Gating const& v) const;
  };
 
 
 
  struct TenTusscher18: public MembraneModelBase<TenTusscher18,18>
  {
    TenTusscher18();
    
    double current(double u, Gating const& v) const;
 
    double current_du(double u, Gating const& v, double h=1e-5) const;
    
    Gating current_dv(double u, Gating const& v) const;
    
    Gating gatingRhs(double u, Gating const& v, double Istim=0) const;
    double gatingRhs1(double u, Gating const& v, int vId, double Istim=0) const;
    
    Gating gatingRhs_du(double u, Gating const& v) const;
    
    GatingJacobian gatingRhs_dv(double u, Gating const& v) const;
  };
 
 
 
  //-----------------------------------------------------------------------------------------------------
 
  template <class Derived, int n, int m>
  struct ActiveStressModelBase
  {
    static int const nState = n;
    
    static int const nTrigger = m;
    
    typedef Dune::FieldVector<double,nState>  State;
    
    typedef Dune::FieldVector<double,nTrigger> Trigger;
    
    typedef Dune::FieldMatrix<double,nState,nState> StateJacobian;
    
    typedef Dune::FieldMatrix<double,nState,nTrigger> StateTriggerJacobian;
    
    ActiveStressModelBase() {}
    
    ActiveStressModelBase(std::string const& name, State const& fix_): nm(name), fix(fix_) {}
    
    std::string const& name() const { return nm; }
    
    State const& restState() const { return fix; }
    
    StateJacobian state_ds(State const& s, Trigger const& t) const
    {
      auto const& me = static_cast<Derived const&>(*this);
      StateJacobian ds;
      State r = me.state(s,t), dr, h;
      for (int i=0; i<nState; ++i)
      {
        h = s;
        h[i] += 1e-5;
        dr = (me.state(h,t)-r)/1e-5;
        for (int j=0; j<nState; ++j)
          ds[j][i] = dr[j];
      }
      return ds;
    }
    
    StateTriggerJacobian state_dt(State const& s, Trigger const& t) const
    {
      auto const& me = static_cast<Derived const&>(*this);
      StateTriggerJacobian ds;
      State r = me.state(s,t), dr;
      Trigger h;
      for (int i=0; i<nTrigger; ++i)
      {
        h = t;
        h[i] += 1e-5;
        dr = (me.state(s,h)-r)/1e-5;
        for (int j=0; j<nState; ++j)
          ds[j][i] = dr[j];
      }
      return ds;
    }
    
    State stress_ds(State const& s) const
    {
      auto const& me = static_cast<Derived const&>(*this);
      State da, h;
      double r = me.stress(s);
      for (int i=0; i<nState; ++i)
      {
        h = s;
        h[i] += 1e-5;
        da[i] = (me.stress(h)-r)/1e-5;
      }
      return da;
    }
    
  private:
    std::string nm;
    State fix;
  };
  
  class XXstress: public ActiveStressModelBase<XXstress,1,1>
  {
  public:
    XXstress():  ActiveStressModelBase<XXstress,1,1>("XXstress",State(0)) {}
    
    double stress(State const& s) const
    {
      return s[0];
    }
    
    State state(State const& s, Trigger const& t) const
    {
      double e0 = .5; 
      double kta = 47.9; 
//       double un = t/120;
      double un = t;
      return e0 * (kta*un-s) / (1+exp(-100*(un-.05))); 
    }
  };
  
  
}
 
#endif