Title here
Summary here
This tutorial was generated from the file projects/PlosOne_mRNA/test/TestSensitivityAnalysisOHaraEndoLiteratePaper.hpp at revision r20146. Note that the code is given in full at the bottom of the page.
PLOS ONE mRNA population study
On this wiki page we describe in detail the code that is used to run the Chaste simulations in the paper
This code, when compiled, runs one of the experimental design files in the test/data folder. This is designed to be done as part of a batch simulation on a supercomputer when generating all results. We recommend running the simulations by first compiling the code using a command such as
scons build=GccOpt compile_only=1 ts=projects/PlosOne_mRNA/test/TestSensitivityAnalysisOHaraEndoLiteratePaper.hpp
We recommend using the GccOptNative or Intel compile options if they are available.
The code is then run using a single thread using, for example:
./${CHASTE}/projects/PlosOne_mRNA/build/optimised/TestSensitivityAnalysisOHaraEndoLiteratePaperRunner –file projects/PlosOne_mRNA/test/data/input_40_surf/exp_design_30_surf_1.dat
this runs and outputs results for all experiments in the file exp_design_30_surf_1.dat
The setting –tag
To generate the results in the paper we ran the simulations with 8 processes per node using the PBS batch scheduler on a cluster to enable execution of multiple files simultaneously.
Figures were generated from the postprocessed output using the matlab files found in the matlab/ folder of this user project.
Code overview
The code only runs if Chaste is set up to use Cvode, hence the #ifdef CHASTE_CVODE line at the top of the file (not shown on wiki, see downloadable project).
Define the header files
#include <cxxtest/TestSuite.h>
#include <ctime>
#include <iostream>
#include <boost/assign/list_of.hpp>
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/lexical_cast.hpp>
#include "ohara_rudy_2011Cvode.hpp"
ohara_rudy_2011Cvode.hpp is autogenerated by chaste_codegen from the file ohara_rudy_2011.cellml in the src/ folder.
#include "AbstractCardiacCell.hpp"
#include "AbstractCvodeCell.hpp"
#include "Exception.hpp"
#include "OutputFileHandler.hpp"
#include "SimpleStimulus.hpp"
#include "EulerIvpOdeSolver.hpp"
#include "RegularStimulus.hpp"
#include "SensitivityDataStructure.hpp"
SensitivityDataStructure.hpp is included from the src/ folder
#include "CellProperties.hpp"
#include "VectorHelperFunctions.hpp"
class TestSensitivityAnalysisOHaraEndoLiteratePaper : public CxxTest::TestSuite
{
public:
/**
* Parameters can be defined at the top of this Test
*/
void TestSensitivityOHaraEndo(void) throw (Exception)
{
Simulates a population of models.
Read in the command line arguments supplied to the executable
//////////// DEFINE PARAMETERS ///////////////
CommandLineArguments* p_args = CommandLineArguments::Instance();
unsigned argc = *(p_args->p_argc); // has the number of arguments.
std::cout << "## " << argc-1 << " arguments supplied.\n" << std::flush;
If no file is set then we display an error explaining valid options
if ( !CommandLineArguments::Instance()->OptionExists("--file"))
{
std::cerr << "TestSensitivityAnalysis::Please input an argument\n"
"* --file the path of an experimental design file\n"
"* --tag tag for identifying output (optional)\n"
"* --store-data record voltage and calcium transients (optional)\n";
return;
}
The FileFinder ensures that the supplied file really does exist. Paths can be either absolute for the system, or relative to the current working directory, which in this case will be the Chaste folder.
FileFinder exp_des_path(CommandLineArguments::Instance()->GetStringCorrespondingToOption("--file"),
RelativeTo::AbsoluteOrCwd);
if (!exp_des_path.Exists())
{
EXCEPTION("No experiment description file found at :" << exp_des_path.GetAbsolutePath());
}
If a tag is supplied then we read it in and append it to the output folder name.
std::string tag = "";
if (CommandLineArguments::Instance()->OptionExists("--tag"))
{
tag=CommandLineArguments::Instance()->GetStringCorrespondingToOption("--tag");
}
std::string foldername = "SensitivityAnalysis/OHara2011_endo" +tag;
Make above directory. This will NOT wipe an existing folder as we have set the false option. If we did not do this then every run of the executable would re-create the folder, deleting any previously calculated results.
OutputFileHandler steady_handler(foldername, false);
Check to see if action potentials and calcium transients are to be stored and set the store_data flag to true if they are.
bool store_data = false;
if (CommandLineArguments::Instance()->OptionExists("--store-data"))
{
store_data=true;
}
The stimulus is set later, we use a simple stimulus for now to set up the simulation. A solver is not required for a Cvode cell so we use a placeholder solver for now.
boost::shared_ptr<SimpleStimulus> p_stimulus(new SimpleStimulus(-80, 0.5, 0));
boost::shared_ptr<EulerIvpOdeSolver> p_solver(new EulerIvpOdeSolver);
Cellohara_rudy_2011FromCellMLCvode* p_model = new Cellohara_rudy_2011FromCellMLCvode(p_solver, p_stimulus);
The following names are fixed and correspond to metadata tags. We record the default parameter values from the model uses and const them to avoid problems!
const double default_g_na = p_model->GetParameter("membrane_fast_sodium_current_conductance");
const double default_g_cal= p_model->GetParameter("membrane_L_type_calcium_current_conductance");
const double default_g_kr = p_model->GetParameter("membrane_rapid_delayed_rectifier_potassium_current_conductance");
const double default_g_ks = p_model->GetParameter("membrane_slow_delayed_rectifier_potassium_current_conductance");
const double default_g_k1 = p_model->GetParameter("membrane_inward_rectifier_potassium_current_conductance");
const double default_j_up = p_model->GetParameter("SR_uptake_current_max");
const double default_g_to = p_model->GetParameter("membrane_transient_outward_current_conductance");
const double default_g_naca = p_model->GetParameter("membrane_sodium_calcium_exchanger_current_conductance");
Increase the max number of steps Cvode can take.
p_model->SetMaxSteps(1e8);
Set the regular stimulus
boost::shared_ptr<RegularStimulus> p_default_stimulus =
boost::static_pointer_cast<RegularStimulus>(p_model->GetStimulusFunction());
Note: s1_period is overwritten later when we enter the stimulation protocol
double s1_period = 1500;
Get the stimulus magnitude from the cellml-derived model
double s_magnitude = p_default_stimulus->GetMagnitude();
Get the stimulus duration from the model
double s_duration = p_default_stimulus->GetDuration();
Time at which the stimulus is applied.
double s_start = 0;
Set the printing time step (ms)
double printing_time_step = 1.0;
Initialise the maximum time step (set later)
double maximum_time_step;
PARAMETERS FOR THE STEADY STATE EXPERIMENTS This is a regular stimulus for steady-state pacing experiments
boost::shared_ptr<RegularStimulus> p_regular_stimulus(new RegularStimulus(s_magnitude, s_duration, s1_period, s_start));
Prevent the solver from missing the stimulus altogether.
if (printing_time_step > s_duration)
{
maximum_time_step = s_duration;
}
else
{
maximum_time_step = printing_time_step;
}
Get the index of the voltage and cytosolic calcium variables in the ODE system so that we can find their solutions for postprocessing later!
boost::shared_ptr<const AbstractOdeSystemInformation> p_ode_info = p_model->GetSystemInformation();
unsigned voltage_index = p_ode_info->GetStateVariableIndex("membrane_voltage");
unsigned calcium_index = p_ode_info->GetStateVariableIndex("cytosolic_calcium_concentration");
Set up a sensitivity data structure (definied in src/SensitivityDataStructure.hpp) to read in all the experiments in the given file.
SensitivityDataStructure scale_factor_data(exp_des_path.GetAbsolutePath());
Create a vector of the pacing cycle lengths we’re interested in.
std::vector<double> pacing_cycle_lengths = boost::assign::list_of(1500)(1000)(900)(800)(700)(600)(500)(450)(400)(350)(300);
/**
* START LOOP OVER EXPERIMENTS
*
* First we print how many experiments are in the file, taken from the SensitivityDataStructure.
*/
std::cout << scale_factor_data.GetNumExperiments() << " experiments in data file.\n";
std::cout << "Running simulations..." << std::endl << std::flush;
Loop over the number of experiments in the file.
for(unsigned experiment_idx = 0; experiment_idx < scale_factor_data.GetNumExperiments(); experiment_idx++)
{
Get the index for the current experiment in the overall experimental design and create the output file name with this index in its name.
unsigned current_experiment = scale_factor_data.GetExperimentNumber(experiment_idx);
std::string current_experiment_string = boost::lexical_cast<std::string>(current_experiment);
std::string OutputFileName;
std::stringstream OutputFileNameStream;
OutputFileNameStream << "DYN_" << current_experiment <<".out";
OutputFileName = OutputFileNameStream.str();
Open the output file
out_stream biomarker_results_file = steady_handler.OpenOutputFile(OutputFileName);
Reset to initial conditions (avoid starting from an unusual state)
p_model->SetStateVariables(p_model->GetInitialConditions());
Initialise the parameters for this experiment.
double gNa_factor = 1;
double gCaL_factor = 1;
double gKr_factor = 1;
double gKs_factor = 1;
double gK1_factor = 1;
double jUp_factor = 1;
double gTo_factor = 1;
double gNaCa_factor = 1;
Parameters are extracted from the SensitivityDataStructure. There is a map between numbers and variable names contained within the function GetScaleFactor in SensitivityDataStructure.cpp. If the parameter is in the file in question then its value is extracted, otherwise not. Note that the parameter names are from the OxMetaData standard used to mark up cellml files. Note that this ‘SensitivityDataStructure’ is different to the one presented in the ‘SimpleSensitivity’ project. If you are looking to do your own analysis, we recommend using the SimpleSensitivty format rather than the one presented here as it more neatly supports parallelisation
if (scale_factor_data.IsParameterInFile("membrane_fast_sodium_current_conductance"))
{
gNa_factor = scale_factor_data.GetScaleFactor(experiment_idx, "membrane_fast_sodium_current_conductance");
}
if (scale_factor_data.IsParameterInFile("membrane_L_type_calcium_current_conductance"))
{
gCaL_factor = scale_factor_data.GetScaleFactor(experiment_idx, "membrane_L_type_calcium_current_conductance");
}
if (scale_factor_data.IsParameterInFile("membrane_rapid_delayed_rectifier_potassium_current_conductance"))
{
gKr_factor = scale_factor_data.GetScaleFactor(experiment_idx, "membrane_rapid_delayed_rectifier_potassium_current_conductance");
}
if (scale_factor_data.IsParameterInFile("membrane_slow_delayed_rectifier_potassium_current_conductance"))
{
gKs_factor = scale_factor_data.GetScaleFactor(experiment_idx, "membrane_slow_delayed_rectifier_potassium_current_conductance");
}
if (scale_factor_data.IsParameterInFile("membrane_inward_rectifier_potassium_current_conductance"))
{
gK1_factor = scale_factor_data.GetScaleFactor(experiment_idx, "membrane_inward_rectifier_potassium_current_conductance");
}
if (scale_factor_data.IsParameterInFile("SR_uptake_current_max"))
{
jUp_factor = scale_factor_data.GetScaleFactor(experiment_idx,"SR_uptake_current_max");
}
if (scale_factor_data.IsParameterInFile("membrane_transient_outward_current_conductance"))
{
gTo_factor = scale_factor_data.GetScaleFactor(experiment_idx,"membrane_transient_outward_current_conductance");
}
if (scale_factor_data.IsParameterInFile("membrane_sodium_calcium_exchanger_current_conductance"))
{
gNaCa_factor = scale_factor_data.GetScaleFactor(experiment_idx,"membrane_sodium_calcium_exchanger_current_conductance");
}
Print the scaling factors to screen
std::cout << "gNa factor = " << gNa_factor << "\n" << std::flush;
std::cout << "gCaL factor = " << gCaL_factor << "\n" << std::flush;
std::cout << "gKr factor = " << gKr_factor << "\n" << std::flush;
std::cout << "gKs factor = " << gKs_factor << "\n" << std::flush;
std::cout << "gK1 factor = " << gK1_factor << "\n" << std::flush;
std::cout << "jUp factor = " << jUp_factor << "\n" << std::flush;
std::cout << "gTo factor = " << gTo_factor << "\n" << std::flush;
std::cout << "gNaCa factor = " << gNaCa_factor << "\n" << std::flush;
The parameter value is set as the default value multiplied by the scaling factor in the file for that parameter.
p_model->SetParameter("membrane_fast_sodium_current_conductance",default_g_na*gNa_factor);
p_model->SetParameter("membrane_L_type_calcium_current_conductance",default_g_cal*gCaL_factor);
p_model->SetParameter("membrane_rapid_delayed_rectifier_potassium_current_conductance",default_g_kr*gKr_factor);
p_model->SetParameter("membrane_slow_delayed_rectifier_potassium_current_conductance",default_g_ks*gKs_factor);
p_model->SetParameter("membrane_inward_rectifier_potassium_current_conductance",default_g_k1*gK1_factor);
p_model->SetParameter("SR_uptake_current_max",default_j_up*jUp_factor);
p_model->SetParameter("membrane_transient_outward_current_conductance",default_g_to*gTo_factor);
p_model->SetParameter("membrane_sodium_calcium_exchanger_current_conductance",default_g_naca*gNaCa_factor);
Use optimisations from #1912 - this causes the cell model variables to be reset.
p_model->SetAutoReset(true);
Start loop over pacing frequencies of interest.
for (unsigned pacing_idx=0; pacing_idx<pacing_cycle_lengths.size(); ++pacing_idx)
{
out_stream data_traces_file;
If the store data flag is set then a file for the data traces is created.
if(store_data)
{
std::stringstream file_name_stream;
file_name_stream << "DataTraces_" << current_experiment << "_" << pacing_cycle_lengths[pacing_idx]<< "ms.dat";
std::string data_traces_file_name = file_name_stream.str();
data_traces_file = steady_handler.OpenOutputFile(data_traces_file_name);
}
/**
* set the pacing period
*/
s1_period = pacing_cycle_lengths[pacing_idx];
p_regular_stimulus->SetPeriod(s1_period);
Assign the regular stimulus to the cell’s stimulus
p_model->SetStimulusFunction(p_regular_stimulus);
Solve the simulation for 98 beats without storing any information
p_model->Solve(0, 98*s1_period, maximum_time_step);
Use optimisations from #1912 - Don’t reset on the next solve call. prevents resetting of variables when solve is called next.
p_model->SetAutoReset(false);
Initialise vectors for biomarkers
std::vector<double> apd80_vector;
std::vector<double> apd30_vector;
std::vector<double> caitd80_vector;
std::vector<double> caitd30_vector;
std::vector<double> calcium_peak_vector;
std::vector<double> delay_vector;
Get the solution for a final 2 paces at this cycle length and get more detail by reducing maximum time step, and print step.
unsigned num_repeats = 2u;
Prepare somewhere to store the traces
std::vector<std::vector<double> > store_results(2*num_repeats+1);
for(unsigned repeats=0; repeats<num_repeats; repeats++)
{
Solve for another beat
OdeSolution solution = p_model->Solve(0, s1_period, maximum_time_step/10, maximum_time_step/10);
Get the solutions for voltage and calcium.
std::vector<double> voltages = solution.GetVariableAtIndex(voltage_index); // Voltage should always be zero
std::vector<double> cytosolic_calcium = solution.GetVariableAtIndex(calcium_index);
If the store_data flag is set to true, save their details into a vector
if (store_data)
{
store_results[0] = solution.rGetTimes(); /* (same on all repeats) */
store_results[2*repeats+1] = voltages;
store_results[2*repeats+2] = cytosolic_calcium;
}
Perform postprocessing on the voltage results
CellProperties voltage_properties(voltages, solution.rGetTimes(),-30);// Threshold for AP is above -70ish.
double apd80, apd30, ap_upstroke_time;//, ap_peak;
try and get the biomarker values of interest
try
{
apd80 = voltage_properties.GetLastActionPotentialDuration(80);
apd30 = voltage_properties.GetLastActionPotentialDuration(30);
ap_upstroke_time = voltage_properties.GetTimeAtLastMaxUpstrokeVelocity();
}
Catch any exceptions thrown
catch (Exception& e)
{
If an AP didn’t occur or repolarisation was not complete then we store the values as -1 to be picked up by the matlab code.
if (e.GetShortMessage()=="AP did not occur, never exceeded threshold voltage." ||
e.GetShortMessage()=="No full action potential was recorded")
{
std::cout << "At Steady pacing of " << s1_period << "ms we did not record any APs\n" << std::flush;
apd80 = -1;
apd30 = -1;
ap_upstroke_time = -1;
}
Otherwise, throw the exception and halt the program
else
{
throw e;
}
}
Calculate a threshold based on max - min of calcium transient for computing the calcium transient durations by looping through the vector
double max_ca = -DBL_MAX;
double min_ca = DBL_MAX;
for (unsigned i=0 ; i<cytosolic_calcium.size(); i++)
{
if (cytosolic_calcium[i] > max_ca)
{
max_ca = cytosolic_calcium[i];
}
if (cytosolic_calcium[i] < min_ca)
{
min_ca = cytosolic_calcium[i];
}
}
double ca_threshold = min_ca + (max_ca - min_ca)*0.2;
Perform postprocessing for the calcium transient
CellProperties calcium_properties(cytosolic_calcium, solution.rGetTimes(),ca_threshold);
Either get the resulting value or flag the failure to trigger as -1 as before.
double caitd80, caitd30, calcium_upstroke_time, calcium_peak;
try
{
caitd80 = calcium_properties.GetLastActionPotentialDuration(80);
caitd30 = calcium_properties.GetLastActionPotentialDuration(30);
calcium_upstroke_time = calcium_properties.GetTimeAtLastMaxUpstrokeVelocity();
calcium_peak = calcium_properties.GetLastPeakPotential();
}
catch (Exception& e)
{
if (e.GetShortMessage()=="AP did not occur, never exceeded threshold voltage." ||
e.GetShortMessage()=="No full action potential was recorded")
{
std::cout << "At Steady pacing of " << s1_period << "ms we did not record any calcium transients\n" << std::flush;
caitd80 = -1;
caitd30 = -1;
calcium_upstroke_time = -1;
calcium_peak = -1;
}
else
{
throw e;
}
}
Calculate the AP-CaT upstroke delay
double delay;
if (calcium_upstroke_time != -1 && ap_upstroke_time!=-1)
{
calcium follows upstroke, so this is positive.
delay = calcium_upstroke_time - ap_upstroke_time;
}
else
{
delay = -1;
}
Store these values into the relevant vector.
apd80_vector.push_back(apd80);
apd30_vector.push_back(apd30);
caitd80_vector.push_back(caitd80);
caitd30_vector.push_back(caitd30);
calcium_peak_vector.push_back(calcium_peak);
delay_vector.push_back(delay);
End loop over the number of repeats for which postprocessing is performed (2 here)
}
Print the voltage and calcium transients to file if required
if(store_data)
{
for (unsigned i=0; i<store_results[0].size(); i++)
{
*data_traces_file << store_results[0][i] << "\t" <<
store_results[1][i] << "\t" <<
store_results[2][i] << "\t" <<
store_results[3][i] << "\t" <<
store_results[4][i] << "\n";
}
data_traces_file->close();
}
Write biomarker results
*biomarker_results_file << pacing_cycle_lengths[pacing_idx] << "\t" <<
apd30_vector[0] << "\t" << apd30_vector[1] << "\t" <<
apd80_vector[0] << "\t" << apd80_vector[1] << "\t" <<
caitd30_vector[0] << "\t" << caitd30_vector[1] << "\t" <<
caitd80_vector[0] << "\t" << caitd80_vector[1] << "\t" <<
delay_vector[0] << "\t" << delay_vector[1] <<"\t" <<
calcium_peak_vector[0] << "\t" << calcium_peak_vector[1] <<"\n" << std::flush;
End loop over pacing frequencies
}
Close the biomarker results file
biomarker_results_file->close();
End the loop over the experiment
}
} /* End of test */
};
#endif //_TESTSENSITIVITYANALYSIS_HPP_
Code
The full code is given below
File name TestSensitivityAnalysisOHaraEndoLiteratePaper.hpp
#include <cxxtest/TestSuite.h>
#include <ctime>
#include <iostream>
#include <boost/assign/list_of.hpp>
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/lexical_cast.hpp>
#include "ohara_rudy_2011Cvode.hpp"
#include "AbstractCardiacCell.hpp"
#include "AbstractCvodeCell.hpp"
#include "Exception.hpp"
#include "OutputFileHandler.hpp"
#include "SimpleStimulus.hpp"
#include "EulerIvpOdeSolver.hpp"
#include "RegularStimulus.hpp"
#include "SensitivityDataStructure.hpp"
#include "CellProperties.hpp"
#include "VectorHelperFunctions.hpp"
class TestSensitivityAnalysisOHaraEndoLiteratePaper : public CxxTest::TestSuite
{
public:
/**
* Parameters can be defined at the top of this Test
*/
void TestSensitivityOHaraEndo(void) throw (Exception)
{
//////////// DEFINE PARAMETERS ///////////////
CommandLineArguments* p_args = CommandLineArguments::Instance();
unsigned argc = *(p_args->p_argc); // has the number of arguments.
std::cout << "## " << argc-1 << " arguments supplied.\n" << std::flush;
if ( !CommandLineArguments::Instance()->OptionExists("--file"))
{
std::cerr << "TestSensitivityAnalysis::Please input an argument\n"
"* --file the path of an experimental design file\n"
"* --tag tag for identifying output (optional)\n"
"* --store-data record voltage and calcium transients (optional)\n";
return;
}
FileFinder exp_des_path(CommandLineArguments::Instance()->GetStringCorrespondingToOption("--file"),
RelativeTo::AbsoluteOrCwd);
if (!exp_des_path.Exists())
{
EXCEPTION("No experiment description file found at :" << exp_des_path.GetAbsolutePath());
}
std::string tag = "";
if (CommandLineArguments::Instance()->OptionExists("--tag"))
{
tag=CommandLineArguments::Instance()->GetStringCorrespondingToOption("--tag");
}
std::string foldername = "SensitivityAnalysis/OHara2011_endo" +tag;
OutputFileHandler steady_handler(foldername, false);
bool store_data = false;
if (CommandLineArguments::Instance()->OptionExists("--store-data"))
{
store_data=true;
}
boost::shared_ptr<SimpleStimulus> p_stimulus(new SimpleStimulus(-80, 0.5, 0));
boost::shared_ptr<EulerIvpOdeSolver> p_solver(new EulerIvpOdeSolver);
Cellohara_rudy_2011FromCellMLCvode* p_model = new Cellohara_rudy_2011FromCellMLCvode(p_solver, p_stimulus);
const double default_g_na = p_model->GetParameter("membrane_fast_sodium_current_conductance");
const double default_g_cal= p_model->GetParameter("membrane_L_type_calcium_current_conductance");
const double default_g_kr = p_model->GetParameter("membrane_rapid_delayed_rectifier_potassium_current_conductance");
const double default_g_ks = p_model->GetParameter("membrane_slow_delayed_rectifier_potassium_current_conductance");
const double default_g_k1 = p_model->GetParameter("membrane_inward_rectifier_potassium_current_conductance");
const double default_j_up = p_model->GetParameter("SR_uptake_current_max");
const double default_g_to = p_model->GetParameter("membrane_transient_outward_current_conductance");
const double default_g_naca = p_model->GetParameter("membrane_sodium_calcium_exchanger_current_conductance");
p_model->SetMaxSteps(1e8);
boost::shared_ptr<RegularStimulus> p_default_stimulus =
boost::static_pointer_cast<RegularStimulus>(p_model->GetStimulusFunction());
double s1_period = 1500;
double s_magnitude = p_default_stimulus->GetMagnitude();
double s_duration = p_default_stimulus->GetDuration();
double s_start = 0;
double printing_time_step = 1.0;
double maximum_time_step;
boost::shared_ptr<RegularStimulus> p_regular_stimulus(new RegularStimulus(s_magnitude, s_duration, s1_period, s_start));
if (printing_time_step > s_duration)
{
maximum_time_step = s_duration;
}
else
{
maximum_time_step = printing_time_step;
}
boost::shared_ptr<const AbstractOdeSystemInformation> p_ode_info = p_model->GetSystemInformation();
unsigned voltage_index = p_ode_info->GetStateVariableIndex("membrane_voltage");
unsigned calcium_index = p_ode_info->GetStateVariableIndex("cytosolic_calcium_concentration");
SensitivityDataStructure scale_factor_data(exp_des_path.GetAbsolutePath());
std::vector<double> pacing_cycle_lengths = boost::assign::list_of(1500)(1000)(900)(800)(700)(600)(500)(450)(400)(350)(300);
/**
* START LOOP OVER EXPERIMENTS
*
* First we print how many experiments are in the file, taken from the SensitivityDataStructure.
*/
std::cout << scale_factor_data.GetNumExperiments() << " experiments in data file.\n";
std::cout << "Running simulations..." << std::endl << std::flush;
for(unsigned experiment_idx = 0; experiment_idx < scale_factor_data.GetNumExperiments(); experiment_idx++)
{
unsigned current_experiment = scale_factor_data.GetExperimentNumber(experiment_idx);
std::string current_experiment_string = boost::lexical_cast<std::string>(current_experiment);
std::string OutputFileName;
std::stringstream OutputFileNameStream;
OutputFileNameStream << "DYN_" << current_experiment <<".out";
OutputFileName = OutputFileNameStream.str();
out_stream biomarker_results_file = steady_handler.OpenOutputFile(OutputFileName);
p_model->SetStateVariables(p_model->GetInitialConditions());
double gNa_factor = 1;
double gCaL_factor = 1;
double gKr_factor = 1;
double gKs_factor = 1;
double gK1_factor = 1;
double jUp_factor = 1;
double gTo_factor = 1;
double gNaCa_factor = 1;
if (scale_factor_data.IsParameterInFile("membrane_fast_sodium_current_conductance"))
{
gNa_factor = scale_factor_data.GetScaleFactor(experiment_idx, "membrane_fast_sodium_current_conductance");
}
if (scale_factor_data.IsParameterInFile("membrane_L_type_calcium_current_conductance"))
{
gCaL_factor = scale_factor_data.GetScaleFactor(experiment_idx, "membrane_L_type_calcium_current_conductance");
}
if (scale_factor_data.IsParameterInFile("membrane_rapid_delayed_rectifier_potassium_current_conductance"))
{
gKr_factor = scale_factor_data.GetScaleFactor(experiment_idx, "membrane_rapid_delayed_rectifier_potassium_current_conductance");
}
if (scale_factor_data.IsParameterInFile("membrane_slow_delayed_rectifier_potassium_current_conductance"))
{
gKs_factor = scale_factor_data.GetScaleFactor(experiment_idx, "membrane_slow_delayed_rectifier_potassium_current_conductance");
}
if (scale_factor_data.IsParameterInFile("membrane_inward_rectifier_potassium_current_conductance"))
{
gK1_factor = scale_factor_data.GetScaleFactor(experiment_idx, "membrane_inward_rectifier_potassium_current_conductance");
}
if (scale_factor_data.IsParameterInFile("SR_uptake_current_max"))
{
jUp_factor = scale_factor_data.GetScaleFactor(experiment_idx,"SR_uptake_current_max");
}
if (scale_factor_data.IsParameterInFile("membrane_transient_outward_current_conductance"))
{
gTo_factor = scale_factor_data.GetScaleFactor(experiment_idx,"membrane_transient_outward_current_conductance");
}
if (scale_factor_data.IsParameterInFile("membrane_sodium_calcium_exchanger_current_conductance"))
{
gNaCa_factor = scale_factor_data.GetScaleFactor(experiment_idx,"membrane_sodium_calcium_exchanger_current_conductance");
}
std::cout << "gNa factor = " << gNa_factor << "\n" << std::flush;
std::cout << "gCaL factor = " << gCaL_factor << "\n" << std::flush;
std::cout << "gKr factor = " << gKr_factor << "\n" << std::flush;
std::cout << "gKs factor = " << gKs_factor << "\n" << std::flush;
std::cout << "gK1 factor = " << gK1_factor << "\n" << std::flush;
std::cout << "jUp factor = " << jUp_factor << "\n" << std::flush;
std::cout << "gTo factor = " << gTo_factor << "\n" << std::flush;
std::cout << "gNaCa factor = " << gNaCa_factor << "\n" << std::flush;
p_model->SetParameter("membrane_fast_sodium_current_conductance",default_g_na*gNa_factor);
p_model->SetParameter("membrane_L_type_calcium_current_conductance",default_g_cal*gCaL_factor);
p_model->SetParameter("membrane_rapid_delayed_rectifier_potassium_current_conductance",default_g_kr*gKr_factor);
p_model->SetParameter("membrane_slow_delayed_rectifier_potassium_current_conductance",default_g_ks*gKs_factor);
p_model->SetParameter("membrane_inward_rectifier_potassium_current_conductance",default_g_k1*gK1_factor);
p_model->SetParameter("SR_uptake_current_max",default_j_up*jUp_factor);
p_model->SetParameter("membrane_transient_outward_current_conductance",default_g_to*gTo_factor);
p_model->SetParameter("membrane_sodium_calcium_exchanger_current_conductance",default_g_naca*gNaCa_factor);
p_model->SetAutoReset(true);
for (unsigned pacing_idx=0; pacing_idx<pacing_cycle_lengths.size(); ++pacing_idx)
{
out_stream data_traces_file;
if(store_data)
{
std::stringstream file_name_stream;
file_name_stream << "DataTraces_" << current_experiment << "_" << pacing_cycle_lengths[pacing_idx]<< "ms.dat";
std::string data_traces_file_name = file_name_stream.str();
data_traces_file = steady_handler.OpenOutputFile(data_traces_file_name);
}
/**
* set the pacing period
*/
s1_period = pacing_cycle_lengths[pacing_idx];
p_regular_stimulus->SetPeriod(s1_period);
p_model->SetStimulusFunction(p_regular_stimulus);
p_model->Solve(0, 98*s1_period, maximum_time_step);
p_model->SetAutoReset(false);
std::vector<double> apd80_vector;
std::vector<double> apd30_vector;
std::vector<double> caitd80_vector;
std::vector<double> caitd30_vector;
std::vector<double> calcium_peak_vector;
std::vector<double> delay_vector;
unsigned num_repeats = 2u;
std::vector<std::vector<double> > store_results(2*num_repeats+1);
for(unsigned repeats=0; repeats<num_repeats; repeats++)
{
OdeSolution solution = p_model->Solve(0, s1_period, maximum_time_step/10, maximum_time_step/10);
std::vector<double> voltages = solution.GetVariableAtIndex(voltage_index); // Voltage should always be zero
std::vector<double> cytosolic_calcium = solution.GetVariableAtIndex(calcium_index);
if (store_data)
{
store_results[0] = solution.rGetTimes(); /* (same on all repeats) */
store_results[2*repeats+1] = voltages;
store_results[2*repeats+2] = cytosolic_calcium;
}
CellProperties voltage_properties(voltages, solution.rGetTimes(),-30);// Threshold for AP is above -70ish.
double apd80, apd30, ap_upstroke_time;//, ap_peak;
try
{
apd80 = voltage_properties.GetLastActionPotentialDuration(80);
apd30 = voltage_properties.GetLastActionPotentialDuration(30);
ap_upstroke_time = voltage_properties.GetTimeAtLastMaxUpstrokeVelocity();
}
catch (Exception& e)
{
if (e.GetShortMessage()=="AP did not occur, never exceeded threshold voltage." ||
e.GetShortMessage()=="No full action potential was recorded")
{
std::cout << "At Steady pacing of " << s1_period << "ms we did not record any APs\n" << std::flush;
apd80 = -1;
apd30 = -1;
ap_upstroke_time = -1;
}
else
{
throw e;
}
}
double max_ca = -DBL_MAX;
double min_ca = DBL_MAX;
for (unsigned i=0 ; i<cytosolic_calcium.size(); i++)
{
if (cytosolic_calcium[i] > max_ca)
{
max_ca = cytosolic_calcium[i];
}
if (cytosolic_calcium[i] < min_ca)
{
min_ca = cytosolic_calcium[i];
}
}
double ca_threshold = min_ca + (max_ca - min_ca)*0.2;
CellProperties calcium_properties(cytosolic_calcium, solution.rGetTimes(),ca_threshold);
double caitd80, caitd30, calcium_upstroke_time, calcium_peak;
try
{
caitd80 = calcium_properties.GetLastActionPotentialDuration(80);
caitd30 = calcium_properties.GetLastActionPotentialDuration(30);
calcium_upstroke_time = calcium_properties.GetTimeAtLastMaxUpstrokeVelocity();
calcium_peak = calcium_properties.GetLastPeakPotential();
}
catch (Exception& e)
{
if (e.GetShortMessage()=="AP did not occur, never exceeded threshold voltage." ||
e.GetShortMessage()=="No full action potential was recorded")
{
std::cout << "At Steady pacing of " << s1_period << "ms we did not record any calcium transients\n" << std::flush;
caitd80 = -1;
caitd30 = -1;
calcium_upstroke_time = -1;
calcium_peak = -1;
}
else
{
throw e;
}
}
double delay;
if (calcium_upstroke_time != -1 && ap_upstroke_time!=-1)
{
delay = calcium_upstroke_time - ap_upstroke_time;
}
else
{
delay = -1;
}
apd80_vector.push_back(apd80);
apd30_vector.push_back(apd30);
caitd80_vector.push_back(caitd80);
caitd30_vector.push_back(caitd30);
calcium_peak_vector.push_back(calcium_peak);
delay_vector.push_back(delay);
}
if(store_data)
{
for (unsigned i=0; i<store_results[0].size(); i++)
{
*data_traces_file << store_results[0][i] << "\t" <<
store_results[1][i] << "\t" <<
store_results[2][i] << "\t" <<
store_results[3][i] << "\t" <<
store_results[4][i] << "\n";
}
data_traces_file->close();
}
*biomarker_results_file << pacing_cycle_lengths[pacing_idx] << "\t" <<
apd30_vector[0] << "\t" << apd30_vector[1] << "\t" <<
apd80_vector[0] << "\t" << apd80_vector[1] << "\t" <<
caitd30_vector[0] << "\t" << caitd30_vector[1] << "\t" <<
caitd80_vector[0] << "\t" << caitd80_vector[1] << "\t" <<
delay_vector[0] << "\t" << delay_vector[1] <<"\t" <<
calcium_peak_vector[0] << "\t" << calcium_peak_vector[1] <<"\n" << std::flush;
}
biomarker_results_file->close();
}
} /* End of test */
};
#endif //_TESTSENSITIVITYANALYSIS_HPP_