Voltage clamping things in Chaste

It's quite useful to be able to clamp the membrane voltage of a cell you're simulating to a particular value. I have found myself wanting to clamp the voltage to a series of values, or even a whole voltage trace from another simulation. There isn't an easy way to do this in Chaste (yet - if I have time I'd like to make it into an easy function), so I'm sharing the bit of code I wrote to make it happen.

If you've not run single-cell cardiac models in Chaste before, look at this tutorial, and if you need to install Chaste on a Debian machine, try this.

First, include ALL THE THINGS:

#include <cxxtest/TestSuite.h>;
#include "ColumnDataReader.hpp"
#include <sstream>
#include <iomanip>
#include "AbstractCardiacCell.hpp"
#include "AbstractCvodeCell.hpp"
#include "CellProperties.hpp"
#include "OutputFileHandler.hpp"
#include "ohara_rudy_2011_endoCvode.hpp"
#include "AbstractParameterisedSystem.hpp"
#include "AbstractCvodeSystem.hpp"
#include <boost /lexical_cast.hpp>
#include "SteadyStateRunner.hpp"
#include "Debug.hpp"
#include "ColumnDataWriter.hpp"

We're clamping a trace to the O'Hara model because that was what I happened to need when I wrote this. This reads in a trace that was output by another Chaste cell simulation in the usual way (see the tutorial linked above), which is stored in projects/BethM/test/TracesForOharaLateSodium/trace.

 

class TestPutTraceIntoOhara : public CxxTest::TestSuite
{
public:
  void TestTakeInTrace() throw(Exception)
   {
    /* Make sure this code is only run if CVODE is installed/enabled on the computer */
#ifdef CHASTE_CVODE

    // open up clamp folder
    std::string folder_root = "projects/BethM/test/TracesForOharaLateSodium/";
                 std::string tracename = "trace";
    FileFinder traces_folder(folder_root + tracename, RelativeTo::ChasteSourceRoot);
                TS_ASSERT_EQUALS(traces_folder.IsDir(), true);
    
    // read in the voltage and time
                boost::shared_ptr<columndatareader> p_reader(new ColumnDataReader(traces_folder,tracename));
                std::vector<double> voltages = p_reader->GetValues("membrane_voltage");
                std::vector</double><double> time_points = p_reader->GetValues("Time");

 

Then as in the tutorial we need to include all the things we need for our simulations: stimulus, solver, time step...

 

    // for running simulations
    double dt = 0.1;

    // stimulus and solver
    double stim_magnitude = -25.5;
    double stim_duration = 3;
    double stim_start = 100;
    double stim_period = 3000;
    boost::shared_ptr<RegularStimulus> p_stimulus(new RegularStimulus(stim_magnitude,stim_duration,stim_period,stim_start));
    boost::shared_ptr<AbstractIvpOdeSolver> p_solver;
    
    // make new model
    boost::shared_ptr<AbstractCvodeCell> p_model(new Cellohara_rudy_2011_endoFromCellMLCvode(p_solver, p_stimulus));
                std::string modelname = p_model->GetSystemName();

 

Because we'll be stopping and starting the simulation at each time step, we won't be able to get the output variables from the solution object, which means we have to make our own data object to store it in.

 

    std::vector<std::vector <double> > outputs;

 

It's usually a good idea to start by running the model to a steady state, just to make sure the variables are what you'd expect - especially if you modify anything.

 

    // run to steady state
    SteadyStateRunner steady_runner(p_model);
    bool result = steady_runner.RunToSteadyState();

 

Then we want to make sure the voltage doesn't do anything by itself between our time points.

 

    // clamp
    p_model->SetVoltageDerivativeToZero();
 

 

Now we run all the way through our input trace and stop and start the model at each time step, changing the voltage each time.

 

          for (unsigned i=0; i < time_points.size()-1; i++)
    {  
      p_model->SetVoltage(voltages[i]);
      double start = time_points[i];
      double end = time_points[i+1];
      p_model->Compute(start,end,dt);

 

Now, we need to get out the time, voltage and any other variables you're interested in.

 

      std::vector<double> state;
      state.push_back(end);
      state.push_back(p_model->GetVoltage());
      state.push_back(p_model->GetAnyVariable("membrane_persistent_sodium_current",end));  
      outputs.push_back(state);
    }

 

What follows is my own messy way of outputting data to files: you might want to use ColumnDataWriter instead, but I don't like the way the directory structure works for that.

 

    
    std::ofstream output_file;
    std::string file_location = getenv("CHASTE_TEST_OUTPUT");
    std::stringstream namestream;
    namestream << file_location <<  "OharaClamp/" << tracename << value <<".dat";
    std::string filename = namestream.str();
    output_file.open(filename.c_str());

    output_file << "Time, membrane voltage, membrane persistent sodium current\n";
    std::vector late_sodiums;
    for (unsigned j= 0; j < outputs.size(); j++)
    {
      output_file << outputs[j][0] << ", " << outputs[j][1] << ", " << outputs[j][2] << "\n";
      late_sodiums.push_back(outputs[j][2]);  
    }
    output_file.close();
  
#else
    /* CVODE is not enable or installed*/
    std::cout << "Cvode is not enabled.\n";
#endif
  }
};

#endif

 

It's a little clunky, as methods go, but it works well enough for my purposes.