Skip to content

Commit

Permalink
Add more sensible threshold for action potentials
Browse files Browse the repository at this point in the history
  • Loading branch information
@mirams
mirams committed Apr 13, 2017
1 parent ea61811 commit fa8a3a4
Show file tree
Hide file tree
Showing 3 changed files with 117 additions and 44 deletions.
42 changes: 41 additions & 1 deletion src/lookup/LookupTableGenerator.cpp
View file
Expand Up @@ -63,6 +63,7 @@ struct ThreadInputData
unsigned mMaxNumPaces;
unsigned mModelIndex;
double mFrequency;
double mVoltageThreshold;
};

template <unsigned DIM>
Expand All @@ -79,7 +80,8 @@ LookupTableGenerator<DIM>::LookupTableGenerator(
mGenerationHasBegun(false),
mMaxRefinementDifference(UNSIGNED_UNSET),
mpParentBox(new ParameterBox<DIM>(NULL)),
mMaxNumPaces(UNSIGNED_UNSET)
mMaxNumPaces(UNSIGNED_UNSET),
mVoltageThreshold(-50.0)
{
// empty
}
Expand Down Expand Up @@ -150,6 +152,11 @@ void LookupTableGenerator<DIM>::GenerateLookupTable()
mUnscaledParameters.push_back(p_model->GetParameter(mParameterNames[i]));
}

// We now do a special run of a model with sodium current set to zero, so we can see the effect
// of simply a stimulus current, and then set the threshold for APs accordingly.
mVoltageThreshold = DetectVoltageThresholdForActionPotential(p_model);
p_model->SetStateVariables(mInitialConditions); // Put the model back to sensible state

// Initial scalings
CornerSet set_of_points = mpParentBox->GetCorners();
assert(set_of_points.size() == pow(2, DIM));
Expand Down Expand Up @@ -263,6 +270,7 @@ void LookupTableGenerator<DIM>::RunEvaluationsForThesePoints(
thread_data[i].mMaxNumPaces = mMaxNumPaces;
thread_data[i].mModelIndex = mModelIndex;
thread_data[i].mFrequency = mFrequency;
thread_data[i].mVoltageThreshold = mVoltageThreshold;

// Launch the ThreadedActionPotential method on this thread
return_code = pthread_create(&threads[i], NULL, ThreadedActionPotential,
Expand Down Expand Up @@ -359,6 +367,7 @@ void* ThreadedActionPotential(void* argument)
ap_runner.SuppressOutput();
ap_runner.SetMaxNumPaces(my_data->mMaxNumPaces);
ap_runner.SetLackOfOneToOneCorrespondenceIsError();
ap_runner.SetVoltageThresholdForRecordingAsActionPotential(my_data->mVoltageThreshold);

// Call the SingleActionPotentialPrediction methods.
try
Expand Down Expand Up @@ -570,6 +579,37 @@ void LookupTableGenerator<DIM>::SetPacingFrequency(double frequency)
mFrequency = frequency;
}

template <unsigned DIM>
double LookupTableGenerator<DIM>::DetectVoltageThresholdForActionPotential(boost::shared_ptr<AbstractCvodeCell> pModel)
{
SingleActionPotentialPrediction ap_runner(pModel);
ap_runner.SuppressOutput();
ap_runner.SetMaxNumPaces(100u);

// We switch off the sodium current and see how high the stimulus makes the voltage go.
if (pModel->HasParameter("membrane_fast_sodium_current_conductance"))
{
const double original_na_conductance = pModel->GetParameter("membrane_fast_sodium_current_conductance");
pModel->SetParameter("membrane_fast_sodium_current_conductance", 0u);

OdeSolution solution = ap_runner.RunSteadyPacingExperiment();

// Put it back where it was! The calling method will reset state variables.
pModel->SetParameter("membrane_fast_sodium_current_conductance", original_na_conductance);

std::vector<double> voltages = solution.GetAnyVariable("membrane_voltage");
double max_voltage = *(std::max_element(voltages.begin(), voltages.end()));
double min_voltage = *(std::min_element(voltages.begin(), voltages.end()));

// Go 10% over the depolarization jump at gNa=0 as a threshold for 'this really is an AP'.
return min_voltage + 1.1 * (max_voltage - min_voltage);
}
else
{
return -50.0; // mV
}
}

#include "SerializationExportWrapperForCpp.hpp"
EXPORT_TEMPLATE_CLASS_SAME_DIMS(LookupTableGenerator)

Expand Down
104 changes: 61 additions & 43 deletions src/lookup/LookupTableGenerator.hpp
View file
Expand Up @@ -36,23 +36,23 @@ OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#ifndef LOOKUPTABLEGENERATOR_HPP_
#define LOOKUPTABLEGENERATOR_HPP_

#include <set>
#include <boost/shared_ptr.hpp>
#include <set>
// Seems that whatever version of ublas I am using now contains boost serialization
// methods for c_vector, which is nice.
#include <boost/serialization/vector.hpp>
#include <boost/serialization/shared_ptr.hpp>
#include <boost/serialization/vector.hpp>
#include "ChasteSerialization.hpp" // Should be included before any other Chaste headers.
#include "ChasteSerializationVersion.hpp"

#include "AbstractCvodeCell.hpp"
#include "ParameterBox.hpp"
#include "Exception.hpp"
#include "UblasVectorInclude.hpp" // Chaste helper header to get c_vectors included with right namespace.
#include "OutputFileHandler.hpp"
#include "SingleActionPotentialPrediction.hpp"
#include "QuantityOfInterest.hpp"
#include "ParameterBox.hpp"
#include "ParameterPointData.hpp"
#include "QuantityOfInterest.hpp"
#include "SingleActionPotentialPrediction.hpp"
#include "UblasVectorInclude.hpp" // Chaste helper header to get c_vectors included with right namespace.

/**
* A class that will generate lookup tables in DIM-dimensional parameter space,
Expand All @@ -68,10 +68,13 @@ OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
* You should add QoIs in order of importance, as the lookup table will be refined for each in turn.
*/
template<unsigned DIM>
template <unsigned DIM>
class LookupTableGenerator
{
private:
private:
/** Needed for testing */
friend class TestLookupTableGenerator;

/** Needed for serialization. */
friend class boost::serialization::access;
/**
Expand All @@ -80,31 +83,35 @@ class LookupTableGenerator
* @param archive the archive
* @param version the current version of this class
*/
template<class Archive>
void serialize(Archive & archive, const unsigned int version)
template <class Archive>
void serialize(Archive& archive, const unsigned int version)
{
archive & mModelIndex;
archive& mModelIndex;
if (version > 0u)
{
archive & mFrequency;
archive& mFrequency;
}
archive& mParameterPoints;
archive& mParameterPointData;
archive& mParameterNames;
archive& mMinimums;
archive& mMaximums;
archive& mUnscaledParameters;
archive& mQuantitiesToRecord;
archive& mQoITolerances;
archive& mMaxNumEvaluations;
archive& mNumEvaluations;
archive& mOutputFileName;
archive& mOutputFolder;
archive& mGenerationHasBegun;
archive& mpParentBox;
archive& mInitialConditions;
archive& mMaxRefinementDifference;
archive& mMaxNumPaces;
if (version > 1u)
{
archive& mVoltageThreshold;
}
archive & mParameterPoints;
archive & mParameterPointData;
archive & mParameterNames;
archive & mMinimums;
archive & mMaximums;
archive & mUnscaledParameters;
archive & mQuantitiesToRecord;
archive & mQoITolerances;
archive & mMaxNumEvaluations;
archive & mNumEvaluations;
archive & mOutputFileName;
archive & mOutputFolder;
archive & mGenerationHasBegun;
archive & mpParentBox;
archive & mInitialConditions;
archive & mMaxRefinementDifference;
archive & mMaxNumPaces;
}

/** Helper wrappers round these long-winded set and iterator names */
Expand Down Expand Up @@ -168,6 +175,9 @@ class LookupTableGenerator
/** The maximum number of paces to do */
unsigned mMaxNumPaces;

/** Threshold to send to action potential evaluation software to say "This is an excited AP" or not. */
double mVoltageThreshold;

/**
* This method will farm out the evaluation of a set of points using multi-threading.
*
Expand All @@ -183,12 +193,21 @@ class LookupTableGenerator
* make sure we can detect and recover the intended state of everything.
*/
LookupTableGenerator()
: mModelIndex(0u),
mpParentBox(NULL)
{};
: mModelIndex(0u),
mpParentBox(NULL){};

public:
/**
* @param pModel a cell model
*
* @return the threshold at which we think a voltage signal is a real
* excited action potential, rather than simply a stimulus current and decay
* so we give the error code NoActionPotential_1 appropriately (rather than really small APDs).
*
* Note this method will mess up state variables and they will need resetting to normal steady state after calling it.
*/
double DetectVoltageThresholdForActionPotential(boost::shared_ptr<AbstractCvodeCell> pModel);

public:
/**
* Constructor
*
Expand Down Expand Up @@ -312,21 +331,20 @@ EXPORT_TEMPLATE_CLASS_SAME_DIMS(LookupTableGenerator)
// Keep track of the archive version we are using
namespace boost
{
namespace serialization
{
/**
namespace serialization
{
/**
* Specify a version number for archive backwards compatibility.
*
* This is how to do BOOST_CLASS_VERSION(AbstractCardiacPde, 1)
* with a templated class.
*/
template <unsigned DIM>
struct version<LookupTableGenerator<DIM> >
{
CHASTE_VERSION_CONTENT(1); // Increment this on serialize method changes.
};
} // namespace serialization
template <unsigned DIM>
struct version<LookupTableGenerator<DIM> >
{
CHASTE_VERSION_CONTENT(2); // Increment this on serialize method changes.
};
} // namespace serialization
} // namespace boost


#endif // LOOKUPTABLEGENERATOR_HPP_
15 changes: 15 additions & 0 deletions test/TestLookupTableGenerator.hpp
View file
Expand Up @@ -164,6 +164,21 @@ class TestLookupTableGenerator : public CxxTest::TestSuite
}
}

void TestVoltageThresholdDetectionAlgorithm() throw(Exception)
{
std::vector<double> thresholds_for_each_model = boost::assign::list_of(-51.3122)(36.8717)(-41.4010)(-41.9142)(-34.6760)(-44.3969);

for (unsigned model_index = 1; model_index < 7; model_index++)
{
SetupModel setup(1.0, model_index); // models at 1 Hz
boost::shared_ptr<AbstractCvodeCell> p_model = setup.GetModel();

LookupTableGenerator<2> generator(model_index, "dummy", "TestLookupTablesThreshold");
double threshold_voltage = generator.DetectVoltageThresholdForActionPotential(p_model);
TS_ASSERT_DELTA(threshold_voltage, thresholds_for_each_model[model_index - 1u], 1e-2);
}
}

void TestLookupTableMaker5d() throw(Exception)
{
unsigned model_index = 2u; // Ten tusscher '06 (table generated for 1 Hz at present)
Expand Down

0 comments on commit fa8a3a4

Please sign in to comment.