Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
HeartConfigRelatedCellFactory.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include "HeartConfigRelatedCellFactory.hpp"
37
38#include <sstream>
39#include "HeartGeometryInformation.hpp"
40#include "ChasteNodesList.hpp"
41#include "HeartFileFinder.hpp"
42#include "CellMLToSharedLibraryConverter.hpp"
43#include "AbstractCardiacCellInterface.hpp"
44#include "Warnings.hpp"
45// This is needed to prevent the chaste_libs=0 build failing
46// on tests that use a dynamically loaded CVODE model
47#include "AbstractCvodeCell.hpp"
48
49template<unsigned SPACE_DIM>
51 : AbstractCardiacCellFactory<SPACE_DIM>(),
52 mDefaultIonicModel(HeartConfig::Instance()->GetDefaultIonicModel())
53{
54 // Read and store possible region definitions
57
58 // Read and store Stimuli
60
61 // if no stimuli provided in XML, need electrodes instead
62 if (mStimuliApplied.size()==0 && (HeartConfig::Instance()->IsElectrodesPresent() == false) )
63 {
64 EXCEPTION("Simulation needs a stimulus (either <Stimuli> or <Electrodes>).");
65 }
66
67 // Read and store Cell Heterogeneities
73
74 // Do we need to convert any CellML files?
76}
77
78template<unsigned SPACE_DIM>
80{
81 if (mDefaultIonicModel.Dynamic().present())
82 {
83 LoadDynamicModel(mDefaultIonicModel, true);
84 }
85 for (unsigned i=0; i<mIonicModelsDefined.size(); i++)
86 {
87 if (mIonicModelsDefined[i].Dynamic().present())
88 {
89 LoadDynamicModel(mIonicModelsDefined[i], true);
90 }
91 }
92}
93
94template<unsigned SPACE_DIM>
96 const cp::ionic_model_selection_type& rModel,
97 bool isCollective)
98{
99 assert(rModel.Dynamic().present());
100 HeartFileFinder file_finder(rModel.Dynamic()->Path());
102 return converter.Convert(file_finder, isCollective);
103}
104
105template<unsigned SPACE_DIM>
109
110template<unsigned SPACE_DIM>
112 boost::shared_ptr<AbstractStimulusFunction> intracellularStimulus,
113 unsigned nodeIndex)
114{
115 cp::ionic_model_selection_type ionic_model = mDefaultIonicModel;
116
117 for (unsigned ionic_model_region_index = 0;
118 ionic_model_region_index < mIonicModelRegions.size();
119 ++ionic_model_region_index)
120 {
121 if (mIonicModelRegions[ionic_model_region_index]->DoesContain(this->GetMesh()->GetNode(nodeIndex)->GetPoint()))
122 {
123 ionic_model = mIonicModelsDefined[ionic_model_region_index];
124 break;
125 }
126 }
127
128 AbstractCardiacCellInterface* p_cell = NULL;
129
130 if (ionic_model.Dynamic().present())
131 {
132#ifndef CHASTE_CAN_CHECKPOINT_DLLS
133 if (HeartConfig::Instance()->GetCheckpointSimulation())
134 {
135 EXCEPTION("Checkpointing is not compatible with dynamically loaded cell models on Mac OS X.");
136 }
137#endif // CHASTE_CAN_CHECKPOINT_DLLS
138 // Load model from shared library
139 DynamicCellModelLoaderPtr p_loader = LoadDynamicModel(ionic_model, false);
140 p_cell = p_loader->CreateCell(this->mpSolver, intracellularStimulus);
141 }
142 else
143 {
144 assert(ionic_model.Hardcoded().present());
145 switch(ionic_model.Hardcoded().get())
146 {
147 case(cp::ionic_models_available_type::LuoRudyI):
148 {
149 p_cell = new CellLuoRudy1991FromCellML(this->mpSolver, intracellularStimulus);
150 break;
151 }
152
153 case(cp::ionic_models_available_type::LuoRudyIBackwardEuler):
154 {
155 p_cell = new CellLuoRudy1991FromCellMLBackwardEuler(this->mpSolver, intracellularStimulus);
156 break;
157 }
158
159 case(cp::ionic_models_available_type::Fox2002):
160 {
161 p_cell = new CellFoxModel2002FromCellML(this->mpSolver, intracellularStimulus);
162 break;
163 }
164
165 case(cp::ionic_models_available_type::Fox2002BackwardEuler):
166 {
167 p_cell = new CellFoxModel2002FromCellMLBackwardEuler(this->mpSolver, intracellularStimulus);
168 break;
169 }
170
171 case(cp::ionic_models_available_type::DifrancescoNoble):
172 {
173 p_cell = new CellDiFrancescoNoble1985FromCellML(this->mpSolver, intracellularStimulus);
174 break;
175 }
176
177 case(cp::ionic_models_available_type::MahajanShiferaw):
178 {
179 p_cell = new CellMahajan2008FromCellML(this->mpSolver, intracellularStimulus);
180 break;
181 }
182
183 case(cp::ionic_models_available_type::MahajanShiferawBackwardEuler):
184 {
185 p_cell = new CellMahajan2008FromCellMLBackwardEuler(this->mpSolver, intracellularStimulus);
186 break;
187 }
188
189 case(cp::ionic_models_available_type::tenTusscher2006):
190 {
191 p_cell = new CellTenTusscher2006EpiFromCellML(this->mpSolver, intracellularStimulus);
192 break;
193 }
194
195 case(cp::ionic_models_available_type::tenTusscher2006BackwardEuler):
196 {
197 p_cell = new CellTenTusscher2006EpiFromCellMLBackwardEuler(this->mpSolver, intracellularStimulus);
198 break;
199 }
200
201 case(cp::ionic_models_available_type::Maleckar):
202 {
203 p_cell = new CellMaleckar2008FromCellML(this->mpSolver, intracellularStimulus);
204 break;
205 }
206
207 case(cp::ionic_models_available_type::HodgkinHuxley):
208 {
209 p_cell = new CellHodgkinHuxley1952FromCellML(this->mpSolver, intracellularStimulus);
210 break;
211 }
212
213 case(cp::ionic_models_available_type::FaberRudy2000):
214 {
215 p_cell = new CellFaberRudy2000FromCellML(this->mpSolver, intracellularStimulus);
216 break;
217 }
218
219 case(cp::ionic_models_available_type::FaberRudy2000Optimised):
220 {
221 p_cell = new CellFaberRudy2000FromCellMLOpt(this->mpSolver, intracellularStimulus);
222 break;
223 }
224
225 default:
226 {
227 //If the ionic model is not in the current enumeration then the XML parser will have picked it up before now!
229 }
230 }
231 }
232
233 // Set parameters
234 try
235 {
236 SetCellParameters(p_cell, nodeIndex);
237 }
238 catch (const Exception& e)
239 {
240 delete p_cell;
241 throw e;
242 }
243 // Generate lookup tables if present
244 p_cell->GetLookupTableCollection();
245
246 return p_cell;
247}
248
249template<unsigned SPACE_DIM>
251 unsigned nodeIndex)
252{
253 // Special case for backwards-compatibility: scale factors
254 for (unsigned ht_index = 0;
255 ht_index < mCellHeterogeneityAreas.size();
256 ++ht_index)
257 {
258 if (mCellHeterogeneityAreas[ht_index]->DoesContain(this->GetMesh()->GetNode(nodeIndex)->GetPoint()))
259 {
260 try
261 {
262 pCell->SetParameter("ScaleFactorGks", mScaleFactorGks[ht_index]);
263 pCell->SetParameter("ScaleFactorGkr", mScaleFactorGkr[ht_index]);
264 pCell->SetParameter("ScaleFactorIto", mScaleFactorIto[ht_index]);
265 }
266 catch (const Exception&)
267 {
268 // Just ignore missing parameter errors in this case
269 }
270 }
271 }
272
274 if (HeartConfig::Instance()->HasDrugDose())
275 {
276 double drug_dose = HeartConfig::Instance()->GetDrugDose();
277 std::map<std::string, std::pair<double, double> > ic50_values = HeartConfig::Instance()->GetIc50Values();
278 for (std::map<std::string, std::pair<double, double> >::iterator it = ic50_values.begin();
279 it != ic50_values.end();
280 ++it)
281 {
282 const std::string param_name = it->first + "_conductance";
283 if (dynamic_cast<AbstractUntemplatedParameterisedSystem*>(pCell)->HasParameter(param_name))
284 {
285 const double original_conductance = pCell->GetParameter(param_name);
286 const double ic50 = it->second.first;
287 const double hill = it->second.second;
288 const double new_conductance = original_conductance/(1.0 + pow(drug_dose/ic50, hill));
289 pCell->SetParameter(param_name, new_conductance);
290 }
291 else
292 {
293 WARNING("Cannot apply drug to cell at node " << nodeIndex << " as it has no parameter named '" << param_name << "'.");
294 }
295 }
296 }
297
298 // SetParameter elements go next so they override the old ScaleFactor* elements.
299 for (unsigned ht_index = 0;
300 ht_index < mCellHeterogeneityAreas.size();
301 ++ht_index)
302 {
303 if (mCellHeterogeneityAreas[ht_index]->DoesContain(this->GetMesh()->GetNode(nodeIndex)->GetPoint()))
304 {
305 for (std::map<std::string, double>::iterator param_it = mParameterSettings[ht_index].begin();
306 param_it != mParameterSettings[ht_index].end();
307 ++param_it)
308 {
309 pCell->SetParameter(param_it->first, param_it->second);
310 }
311 }
312 }
313}
314
315template<unsigned SPACE_DIM>
317 unsigned nodeIndex)
318{
319 boost::shared_ptr<MultiStimulus> node_specific_stimulus(new MultiStimulus());
320 // Check which of the defined stimuli contain the current node
321 for (unsigned stimulus_index = 0;
322 stimulus_index < mStimuliApplied.size();
323 ++stimulus_index)
324 {
325 if (mStimulatedAreas[stimulus_index]->DoesContain(this->GetMesh()->GetNode(nodeIndex)->GetPoint()))
326 {
327 node_specific_stimulus->AddStimulus(mStimuliApplied[stimulus_index]);
328 }
329 }
330 pCell->SetIntracellularStimulusFunction(node_specific_stimulus);
331}
332
333template<unsigned SPACE_DIM>
335{
336 boost::shared_ptr<MultiStimulus> node_specific_stimulus(new MultiStimulus());
337
338 // Check which of the defined stimuli contain the current node
339 for (unsigned stimulus_index = 0;
340 stimulus_index < mStimuliApplied.size();
341 ++stimulus_index)
342 {
343 if (mStimulatedAreas[stimulus_index]->DoesContain(pNode->GetPoint()))
344 {
345 node_specific_stimulus->AddStimulus(mStimuliApplied[stimulus_index]);
346 }
347 }
348
349 unsigned node_index = pNode->GetIndex();
350 return CreateCellWithIntracellularStimulus(node_specific_stimulus, node_index);
351}
352
353// LCOV_EXCL_START
354template<unsigned SPACE_DIM>
359// LCOV_EXCL_STOP
360
361template<>
363{
364 std::string mesh_file_name = HeartConfig::Instance()->GetMeshName();
365 //files containing list of nodes on each surface
366 std::string epi_surface = mesh_file_name + ".epi";
367 std::string lv_surface = mesh_file_name + ".lv";
368 std::string rv_surface = mesh_file_name + ".rv";
369
370
371 //create the HeartGeometryInformation object
372 //HeartGeometryInformation<3u> info(mesh, epi_surface, lv_surface, rv_surface, true);
373 HeartGeometryInformation<3u> info(*(this->GetMesh()), epi_surface, lv_surface, rv_surface, true);
374
375 //We need the fractions of epi and endo layer supplied by the user
376 double epi_fraction = HeartConfig::Instance()->GetEpiLayerFraction();
377 double endo_fraction = HeartConfig::Instance()->GetEndoLayerFraction();
378
379 //given the fraction of each layer, compute the distance map and fill in the vector
380 info.DetermineLayerForEachNode(epi_fraction,endo_fraction);
381 //get the big heterogeneity vector
382 std::vector<unsigned> heterogeneity_node_list;
383 for (unsigned index=0; index<this->GetMesh()->GetNumNodes(); index++)
384 {
385 heterogeneity_node_list.push_back(info.rGetLayerForEachNode()[index]);
386 }
387
388 std::vector<Node<3u>*> epi_nodes;
389 std::vector<Node<3u>*> mid_nodes;
390 std::vector<Node<3u>*> endo_nodes;
391
392 //create the list of (pointer to object) nodes in each layer from the heterogeneities vector that was just filled in
393 for (unsigned node_index = 0; node_index < this->GetMesh()->GetNumNodes(); node_index++)
394 {
395 if (this->GetMesh()->GetDistributedVectorFactory()->IsGlobalIndexLocal(node_index) )
396 {
397 switch (heterogeneity_node_list[node_index])
398 {
399 //epi
400 case 2u:
401 {
402 epi_nodes.push_back(this->GetMesh()->GetNode(node_index));
403 break;
404 }
405 //mid
406 case 1u:
407 {
408 mid_nodes.push_back(this->GetMesh()->GetNode(node_index));
409 break;
410 }
411 //endo
412 case 0u:
413 {
414 endo_nodes.push_back(this->GetMesh()->GetNode(node_index));
415 break;
416 }
417 default:
419 }
420 }
421 }
422 //assert((endo_nodes.size()+epi_nodes.size()+mid_nodes.size())==this->GetMesh()->GetNumNodes());
423
424 // now the 3 list of pointer to nodes need to be pushed into the mCellHeterogeneityAreas vector,
425 // IN THE ORDER PRESCRIBED BY THE USER IN THE XML FILE!
426 // This is because the corresponding scale factors are already read in that order.
427
428 //these three unsigned tell us in which order the user supplied each layer in the XML file
429 unsigned user_supplied_epi_index = HeartConfig::Instance()->GetEpiLayerIndex();
430 unsigned user_supplied_mid_index = HeartConfig::Instance()->GetMidLayerIndex();
431 unsigned user_supplied_endo_index = HeartConfig::Instance()->GetEndoLayerIndex();
432
433 //these three should have been set to 0, 1 and 2 by HeartConfig::GetCellHeterogeneities
434 assert(user_supplied_epi_index<3);
435 assert(user_supplied_mid_index<3);
436 assert(user_supplied_endo_index<3);
437
438 //pute them in a vector
439 std::vector<unsigned> user_supplied_indices;
440 user_supplied_indices.push_back(user_supplied_epi_index);
441 user_supplied_indices.push_back(user_supplied_mid_index);
442 user_supplied_indices.push_back(user_supplied_endo_index);
443
444 //figure out who goes first
445
446 //loop three times
447 for (unsigned layer_index=0; layer_index<3; layer_index++)
448 {
449 unsigned counter = 0;
450 //find the corresponding index
451 for (unsigned supplied_index = 0; supplied_index<user_supplied_indices.size(); supplied_index++)
452 {
453 if (user_supplied_indices[supplied_index] == layer_index)
454 {
455 break;
456 }
457 counter++;
458 }
459
460 //create the node lists based on the calculations above
461 if (counter==0)
462 {
463 mCellHeterogeneityAreas.push_back(boost::shared_ptr<AbstractChasteRegion<3u> >(new ChasteNodesList<3u>(epi_nodes)) );
464 }
465 if (counter==1)
466 {
467 mCellHeterogeneityAreas.push_back(boost::shared_ptr<AbstractChasteRegion<3u> >(new ChasteNodesList<3u>(mid_nodes)) );
468 }
469 if (counter==2)
470 {
471 mCellHeterogeneityAreas.push_back(boost::shared_ptr<AbstractChasteRegion<3u> >(new ChasteNodesList<3u>(endo_nodes)) );
472 }
473 assert(counter<3);
474 }
475 assert(mCellHeterogeneityAreas.size()==3);
476}
477
478// Explicit instantiation
#define EXCEPTION(message)
#define NEVER_REACHED
void SetIntracellularStimulusFunction(boost::shared_ptr< AbstractStimulusFunction > pStimulus)
virtual void SetParameter(const std::string &rParameterName, double value)=0
virtual double GetParameter(const std::string &rParameterName)=0
virtual AbstractLookupTableCollection * GetLookupTableCollection()
DynamicCellModelLoaderPtr Convert(const FileFinder &rFilePath, bool isCollective=true)
std::vector< cp::ionic_model_selection_type > mIonicModelsDefined
DynamicCellModelLoaderPtr LoadDynamicModel(const cp::ionic_model_selection_type &rModel, bool isCollective)
std::vector< boost::shared_ptr< AbstractChasteRegion< SPACE_DIM > > > mIonicModelRegions
std::vector< std::map< std::string, double > > mParameterSettings
AbstractCardiacCellInterface * CreateCardiacCellForTissueNode(Node< SPACE_DIM > *pNode)
void SetCellIntracellularStimulus(AbstractCardiacCellInterface *pCell, unsigned nodeIndex)
std::vector< boost::shared_ptr< AbstractChasteRegion< SPACE_DIM > > > mCellHeterogeneityAreas
AbstractCardiacCellInterface * CreateCellWithIntracellularStimulus(boost::shared_ptr< AbstractStimulusFunction > intracellularStimulus, unsigned nodeIndex)
std::vector< boost::shared_ptr< AbstractChasteRegion< SPACE_DIM > > > mStimulatedAreas
std::vector< boost::shared_ptr< AbstractStimulusFunction > > mStimuliApplied
void SetCellParameters(AbstractCardiacCellInterface *pCell, unsigned nodeIndex)
void GetStimuli(std::vector< boost::shared_ptr< AbstractStimulusFunction > > &rStimuliApplied, std::vector< boost::shared_ptr< AbstractChasteRegion< DIM > > > &rStimulatedAreas) const
std::string GetMeshName() const
double GetDrugDose() const
bool IsElectrodesPresent() const
unsigned GetMidLayerIndex()
void GetIonicModelRegions(std::vector< boost::shared_ptr< AbstractChasteRegion< DIM > > > &rDefinedRegions, std::vector< cp::ionic_model_selection_type > &rIonicModels) const
unsigned GetEpiLayerIndex()
double GetEndoLayerFraction()
std::map< std::string, std::pair< double, double > > GetIc50Values()
double GetEpiLayerFraction()
unsigned GetEndoLayerIndex()
void GetCellHeterogeneities(std::vector< boost::shared_ptr< AbstractChasteRegion< DIM > > > &rCellHeterogeneityRegions, std::vector< double > &rScaleFactorGks, std::vector< double > &rScaleFactorIto, std::vector< double > &rScaleFactorGkr, std::vector< std::map< std::string, double > > *pParameterSettings)
static HeartConfig * Instance()
Definition Node.hpp:59
ChastePoint< SPACE_DIM > GetPoint() const
Definition Node.cpp:133
unsigned GetIndex() const
Definition Node.cpp:158