#include <NhsModelWithBackwardSolver.hpp>

Inheritance diagram for NhsModelWithBackwardSolver:

Collaboration diagram for NhsModelWithBackwardSolver:

Public Member Functions
	NhsModelWithBackwardSolver ()

void	RunDoNotUpdate (double startTime, double endTime, double timestep)

double	GetNextActiveTension ()

void	RunAndUpdate (double startTime, double endTime, double timestep)

Public Member Functions inherited from NhsContractionModel
	NhsContractionModel ()

void	SetStretchAndStretchRate (double lambda, double dlambdaDt)

void	SetInputParameters (ContractionModelInputParameters &rInputParameters)

void	SetIntracellularCalciumConcentration (double calciumConcentration)

double	GetCalciumTroponinValue ()

void	EvaluateYDerivatives (double time, const std::vector< double > &rY, std::vector< double > &rDY)

double	GetActiveTension ()

bool	IsStretchDependent ()

bool	IsStretchRateDependent ()

Public Member Functions inherited from AbstractOdeBasedContractionModel
	AbstractOdeBasedContractionModel (unsigned numStateVariables)

void	UpdateStateVariables ()

Public Member Functions inherited from AbstractOdeSystem
	AbstractOdeSystem (unsigned numberOfStateVariables)

virtual	~AbstractOdeSystem ()

virtual bool	CalculateStoppingEvent (double time, const std::vector< double > &rY)

virtual double	CalculateRootFunction (double time, const std::vector< double > &rY)

bool	GetUseAnalyticJacobian ()

const std::vector< double > &	rGetConstStateVariables () const

Public Member Functions inherited from AbstractParameterisedSystem< std::vector< double > >
	AbstractParameterisedSystem (unsigned numberOfStateVariables)

std::vector< double > &	rGetStateVariables ()

std::vector< double >	GetStateVariables ()

void	SetStateVariables (const std::vector< double > &rStateVariables)

double	GetStateVariable (unsigned index) const

double	GetStateVariable (const std::string &rName) const

void	SetStateVariable (unsigned index, double newValue)

void	SetStateVariable (const std::string &rName, double newValue)

virtual void	VerifyStateVariables ()

void	SetDefaultInitialConditions (const std::vector< double > &rInitialConditions)

void	SetDefaultInitialCondition (unsigned index, double initialCondition)

std::vector< double >	GetInitialConditions () const

void	ResetToInitialConditions ()

double	GetParameter (unsigned index) const

double	GetParameter (const std::string &rName) const

void	SetParameter (const std::string &rName, double value)

void	SetParameter (unsigned index, double value)

double	GetAnyVariable (unsigned index, double time=0.0, std::vector< double > *pDerivedQuantities=NULL)

double	GetAnyVariable (const std::string &rName, double time=0.0, std::vector< double > *pDerivedQuantities=NULL)

void	SetAnyVariable (unsigned index, double value)

void	SetAnyVariable (const std::string &rName, double value)

virtual std::vector< double >	ComputeDerivedQuantities (double time, const std::vector< double > &rState)

std::vector< double >	ComputeDerivedQuantitiesFromCurrentState (double time)

Public Member Functions inherited from AbstractUntemplatedParameterisedSystem
	AbstractUntemplatedParameterisedSystem (unsigned numberOfStateVariables)

virtual	~AbstractUntemplatedParameterisedSystem ()

boost::shared_ptr< const AbstractOdeSystemInformation >	GetSystemInformation () const

std::string	GetSystemName () const

unsigned	GetNumberOfAttributes () const

bool	HasAttribute (const std::string &rName) const

double	GetAttribute (const std::string &rName) const

unsigned	GetNumberOfStateVariables () const

const std::vector< std::string > &	rGetStateVariableNames () const

const std::vector< std::string > &	rGetStateVariableUnits () const

unsigned	GetStateVariableIndex (const std::string &rName) const

bool	HasStateVariable (const std::string &rName) const

std::string	GetStateVariableUnits (unsigned index) const

unsigned	GetNumberOfParameters () const

const std::vector< std::string > &	rGetParameterNames () const

const std::vector< std::string > &	rGetParameterUnits () const

unsigned	GetParameterIndex (const std::string &rName) const

bool	HasParameter (const std::string &rName) const

std::string	GetParameterUnits (unsigned index) const

unsigned	GetNumberOfDerivedQuantities () const

const std::vector< std::string > &	rGetDerivedQuantityNames () const

const std::vector< std::string > &	rGetDerivedQuantityUnits () const

unsigned	GetDerivedQuantityIndex (const std::string &rName) const

bool	HasDerivedQuantity (const std::string &rName) const

std::string	GetDerivedQuantityUnits (unsigned index) const

unsigned	GetAnyVariableIndex (const std::string &rName) const

bool	HasAnyVariable (const std::string &rName) const

std::string	GetAnyVariableUnits (unsigned index) const

std::string	GetAnyVariableUnits (const std::string &rName) const

Public Member Functions inherited from AbstractContractionModel
	AbstractContractionModel ()

virtual	~AbstractContractionModel ()

void	SetStretch (double stretch)

Private Member Functions
double	ImplicitSolveForQ ()

void	CalculateCaTropAndZDerivatives (double calciumTroponin, double z, double Q, double &dCaTrop, double &dz)

void	CalculateBackwardEulerResidual (double calciumTroponin, double z, double Q, double &residualComponent1, double &residualComponent2)

Private Attributes
double	mDt

Static Private Attributes
static const double	mTolerance = 1e-10

Additional Inherited Members
Protected Member Functions inherited from NhsContractionModel
void	CalculateCalciumTrop50 ()

double	CalculateT0 (double z)

Protected Member Functions inherited from AbstractParameterisedSystem< std::vector< double > >
std::string	DumpState (const std::string &rMessage)

std::string	DumpState (const std::string &rMessage, std::vector< double >Y)

std::string	DumpState (const std::string &rMessage, std::vector< double >Y, double time)

void	CheckParametersOnLoad (const std::vector< double > &rParameters, const std::vector< std::string > &rParameterNames)

Protected Attributes inherited from NhsContractionModel
double	mLambda

double	mDLambdaDt

double	mCalciumI

double	mCalciumTrop50

double	mK1

double	mK2

Protected Attributes inherited from AbstractOdeBasedContractionModel
std::vector< double >	mTemporaryStateVariables

double	mTime

Protected Attributes inherited from AbstractOdeSystem
bool	mUseAnalyticJacobian

Protected Attributes inherited from AbstractParameterisedSystem< std::vector< double > >
std::vector< double >	mStateVariables

std::vector< double >	mParameters

Protected Attributes inherited from AbstractUntemplatedParameterisedSystem
unsigned	mNumberOfStateVariables

boost::shared_ptr < AbstractOdeSystemInformation >	mpSystemInfo

Static Protected Attributes inherited from NhsContractionModel
static const double	mKon = 100

static const double	mKrefoff = 0.2

static const double	mGamma = 2

static const double	mCalciumTroponinMax = 0.07

static const double	mAlphaR1 = 0.002

static const double	mAlphaR2 = 0.0017

static const double	mKZ = 0.15

static const unsigned	mNr = 3u

static const double	mBeta1 = -4

static const double	mAlpha0 = 0.008

static const unsigned	mN = 3u

static const double	mZp = 0.85

static const double	mCalcium50ref = 0.00105

static const double	mTref = 56.2

static const double	mBeta0 = 4.9

static const double	mA = 0.35

static const double	mA1 = -29

static const double	mA2 = 138

static const double	mA3 = 129

static const double	mAlpha1 = 0.03

static const double	mAlpha2 = 0.130

static const double	mAlpha3 = 0.625

Detailed Description

The full backward Euler method on the NHS system compute the solution w^{n+1} in R^5 by: w^{n+1} - dt g(w^{n+1}) = w^n where the ODE system is dw/dt = g(w).

This would involve solving a 5d nonlinear system at each timestep. However, only g1 and g2 (ie dCatrop/dt and dz/dt) are nonlinear, g3, g4 and g5 (corresponding to dQi/dt) are linear and uncoupled. Therefore the backward euler solutions of Q1,Q2,Q3 can be computed immediately, leaving a 2D nonlinear system to be solved using newton's method

Definition at line 56 of file NhsModelWithBackwardSolver.hpp.

Constructor & Destructor Documentation

NhsModelWithBackwardSolver::NhsModelWithBackwardSolver ( )

Constructor

Definition at line 106 of file NhsModelWithBackwardSolver.cpp.

References AbstractOdeBasedContractionModel::mTemporaryStateVariables.

Member Function Documentation

void NhsModelWithBackwardSolver::CalculateBackwardEulerResidual	(	double	calciumTroponin,
		double	z,
		double	Q,
		double &	residualComponent1,
		double &	residualComponent2
	)

private

Compute the residual function for the 2D nonlinear system when Backward Euler is used The backward Euler discretisation is: w^{n+1} - dt g(w^{n+1}) = w^n where: w = (ca_trop, z) and the ODE system is: dw/dt = g(w)

so the residual is: f = w^{n+1} - dt g(w^{n+1}) - w^n

Parameters

calciumTroponin	Current guess for Ca_trop value
z	current guess for z
Q	= Q1+Q2+Q3, where Qi already computed at next timestep
residualComponent1	Returned value - first component of residual
residualComponent2	Returned value - second component of residual

Definition at line 95 of file NhsModelWithBackwardSolver.cpp.

References CalculateCaTropAndZDerivatives(), mDt, and AbstractOdeBasedContractionModel::mTemporaryStateVariables.

Referenced by RunDoNotUpdate().

void NhsModelWithBackwardSolver::CalculateCaTropAndZDerivatives	(	double	calciumTroponin,
		double	z,
		double	Q,
		double &	dCaTrop,
		double &	dz
	)

private

The same as EvaluateYDerivatives in NhsContractionModel, but doesn't use std::vectors (for efficiency, and because this class is hardcoded and hand-optimised for backward euler), and just returns the derivatives of the first two components

Parameters

calciumTroponin	Calcium troponin
z	z
Q	Q=Q1+Q2+Q3
dCaTrop	the returned value of dCaTrop/dt
dz	the returned value of dz/dt

Definition at line 54 of file NhsModelWithBackwardSolver.cpp.

References NhsContractionModel::CalculateT0(), EXCEPTION, NhsContractionModel::mA, NhsContractionModel::mAlpha0, NhsContractionModel::mAlphaR1, NhsContractionModel::mAlphaR2, NhsContractionModel::mCalciumI, NhsContractionModel::mCalciumTrop50, NhsContractionModel::mCalciumTroponinMax, NhsContractionModel::mGamma, NhsContractionModel::mKon, NhsContractionModel::mKrefoff, NhsContractionModel::mKZ, NhsContractionModel::mN, NhsContractionModel::mNr, and NhsContractionModel::mTref.

Referenced by CalculateBackwardEulerResidual().

double NhsModelWithBackwardSolver::GetNextActiveTension ( )

virtual

Returns: the active tension corresponding to the stored state variables computed from the last RunDoNotUpdate(), ie the active tension at the next time. Note that calling GetActiveTension() on the base class will use the internal state variables and return the active tension at the last time, if RunDoNotUpdate() has been called but UpdateStateVariables() has not

Reimplemented from NhsContractionModel.

Definition at line 180 of file NhsModelWithBackwardSolver.cpp.

References NhsContractionModel::CalculateT0(), NhsContractionModel::mA, and AbstractOdeBasedContractionModel::mTemporaryStateVariables.

double NhsModelWithBackwardSolver::ImplicitSolveForQ ( )

private

Solve for Q1,Q2,Q3 (and therefore Q) implicitly using backward euler. These can be done directly as the rhs is linear in Qi

Returns: Q=Q1+Q2+Q3

Definition at line 45 of file NhsModelWithBackwardSolver.cpp.

References NhsContractionModel::mA1, NhsContractionModel::mA2, NhsContractionModel::mA3, NhsContractionModel::mAlpha1, NhsContractionModel::mAlpha2, NhsContractionModel::mAlpha3, NhsContractionModel::mDLambdaDt, mDt, and AbstractOdeBasedContractionModel::mTemporaryStateVariables.

Referenced by RunDoNotUpdate().

void NhsModelWithBackwardSolver::RunAndUpdate	(	double	startTime,
		double	endTime,
		double	timestep
	)

virtual

Overload the RunAndUpdate() method too, as that would use the base class's default (euler) solver

Parameters

startTime
endTime
timestep

Reimplemented from AbstractOdeBasedContractionModel.

Definition at line 195 of file NhsModelWithBackwardSolver.cpp.

References RunDoNotUpdate(), and AbstractOdeBasedContractionModel::UpdateStateVariables().

void NhsModelWithBackwardSolver::RunDoNotUpdate	(	double	startTime,
		double	endTime,
		double	timestep
	)

virtual

Solves for the new state variables at the given end time using the implicit method. Note that the internal state variables are not altered, the solution is saved instead. Call UpdateStateVariables() to update, and GetNextActiveTension() to get the solved active tension

The state variables are not updated because this solve will be called as part of the newton iteration (ie guess stretch, see what the new active tension is) in a fully implicit method.

Note: overloaded from the method in AbstractOdeBasedContractionModel, which just does a simple Euler solve

Parameters