Kerchoffs2003ContractionModel.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 
00030 #include "Kerchoffs2003ContractionModel.hpp"
00031 #include "Exception.hpp"
00032 #include "TimeStepper.hpp"
00033 #include <iostream>
00034 
00035 const double Kerchoffs2003ContractionModel::a6 = 2.0; // 1/um
00036 const double Kerchoffs2003ContractionModel::a7 = 1.5; // um
00037 const double Kerchoffs2003ContractionModel::T0 = 180; // kPa
00038 const double Kerchoffs2003ContractionModel::Ea = 20;  // 1/um
00039 const double Kerchoffs2003ContractionModel::v0 = 0.0075; // um/ms
00040 const double Kerchoffs2003ContractionModel::ls0 = 1.9; // um
00041 const double Kerchoffs2003ContractionModel::tr = 75; // ms
00042 const double Kerchoffs2003ContractionModel::td = 75; // ms
00043 const double Kerchoffs2003ContractionModel::b = 150; // ms/um
00044 const double Kerchoffs2003ContractionModel::ld = -0.4; // um
00045 
00046 Kerchoffs2003ContractionModel::Kerchoffs2003ContractionModel()
00047     : AbstractOdeBasedContractionModel(1)
00048 {
00049     mpSystemInfo = OdeSystemInformation<Kerchoffs2003ContractionModel>::Instance();
00050 
00051     mSarcomereLength = ls0;
00052 
00053     mStateVariables.push_back(mSarcomereLength-1.0/Ea); //steady state
00054 
00055     mIsActivated = false;
00056     mElectricallyUnactivated = true;
00057     mActivationTime = 0.0;
00058     mTime = 0.0;
00059 }
00060 
00061 
00062 void Kerchoffs2003ContractionModel::EvaluateYDerivatives(double time,
00063                                                          const std::vector<double>& rY,
00064                                                          std::vector<double>& rDY)
00065 {
00066     double lc = rY[0];
00067     rDY[0]=( Ea*(mSarcomereLength-lc) - 1 )*v0;
00068 }
00069 
00070 
00071 void Kerchoffs2003ContractionModel::SetInputParameters(ContractionModelInputParameters& rInputParameters)
00072 {
00073     assert(rInputParameters.voltage != DOUBLE_UNSET);
00074 
00075     if (mIsActivated && (rInputParameters.voltage < mDeactivationVoltage))
00076     {
00077         // inactive (resting) - note don't set mIsActivated=false yet 
00078         // as the cell may yet be producing force, and the code is such
00079         // that if mIsActivated=false, Ta=0
00080         mElectricallyUnactivated = true;
00081     }
00082 
00083     if (!mIsActivated && (rInputParameters.voltage > mActivationVoltage))
00084     {
00085         // activated
00086         mIsActivated = true;
00087         mActivationTime = mTime;
00088         mElectricallyUnactivated = false;
00089     }
00090 }
00091 
00092 void Kerchoffs2003ContractionModel::SetStretchAndStretchRate(double stretch, double stretchRate)
00093 {
00094     mSarcomereLength = stretch*ls0;
00095 }
00096 
00097 
00098 double Kerchoffs2003ContractionModel::GetActiveTension(double lc)
00099 {
00100     double f_iso = 0;
00101     if(lc > a7)
00102     {
00103         f_iso = T0 * pow(tanh(a6*(lc-a7)),2);
00104     }
00105 
00106     double f_twitch = 0;
00107     double t_max = b*(mSarcomereLength - ld);
00108     if(mIsActivated)
00109     {
00110         double t_a = mTime - mActivationTime;
00111 
00112         if(t_a < t_max)
00113         {
00114             f_twitch = pow( tanh(t_a/tr)*tanh((t_max-t_a)/td), 2);
00115         }
00116         else if(mElectricallyUnactivated)
00117         {
00118             // t_a < t_ma => f_twitch=0 => Ta=0
00119             // In this case, if electrically unactivated as well, 
00120             // set the state to be unactivated.
00121             mIsActivated = false;
00122         }
00123     }
00124 
00125     // expl is unstable for dt = 0.01, 0.001, impl is fine
00126     return (mSarcomereLength/ls0)*f_iso*f_twitch*(mSarcomereLength-lc)*Ea;
00127 }
00128 
00129 
00130 double Kerchoffs2003ContractionModel::GetActiveTension()
00131 {
00132     return GetActiveTension(mStateVariables[0]);
00133 }
00134 
00135 double Kerchoffs2003ContractionModel::GetNextActiveTension()
00136 {
00137     return GetActiveTension(mTemporaryStateVariables[0]);
00138 }
00139 
00140 template<>
00141 void OdeSystemInformation<Kerchoffs2003ContractionModel>::Initialise()
00142 {
00143     this->mVariableNames.push_back("lc");
00144     this->mVariableUnits.push_back("um");
00145     this->mInitialised = true;
00146 }

Generated on Tue May 31 14:31:41 2011 for Chaste by  doxygen 1.5.5