Chaste  Release::2018.1
HeartConfigRelatedCellFactory.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, 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 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF 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 
49 template<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 
78 template<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 
94 template<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 
105 template<unsigned SPACE_DIM>
107 {
108 }
109 
110 template<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
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 
249 template<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 
315 template<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 
333 template<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
354 template<unsigned SPACE_DIM>
356 {
358 }
359 // LCOV_EXCL_STOP
360 
361 template<>
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
479 template class HeartConfigRelatedCellFactory<1u>;
480 template class HeartConfigRelatedCellFactory<2u>;
481 template class HeartConfigRelatedCellFactory<3u>;
const std::vector< HeartLayerType > & rGetLayerForEachNode()
void DetermineLayerForEachNode(double epiFraction, double endoFraction)
void GetIonicModelRegions(std::vector< boost::shared_ptr< AbstractChasteRegion< DIM > > > &rDefinedRegions, std::vector< cp::ionic_model_selection_type > &rIonicModels) const
Definition: Node.hpp:58
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)
virtual void SetParameter(const std::string &rParameterName, double value)=0
DynamicCellModelLoaderPtr LoadDynamicModel(const cp::ionic_model_selection_type &rModel, bool isCollective)
#define EXCEPTION(message)
Definition: Exception.hpp:143
std::vector< boost::shared_ptr< AbstractChasteRegion< SPACE_DIM > > > mCellHeterogeneityAreas
unsigned GetEndoLayerIndex()
std::vector< boost::shared_ptr< AbstractStimulusFunction > > mStimuliApplied
void SetIntracellularStimulusFunction(boost::shared_ptr< AbstractStimulusFunction > pStimulus)
AbstractCardiacCellInterface * CreateCellWithIntracellularStimulus(boost::shared_ptr< AbstractStimulusFunction > intracellularStimulus, unsigned nodeIndex)
#define NEVER_REACHED
Definition: Exception.hpp:206
std::vector< boost::shared_ptr< AbstractChasteRegion< SPACE_DIM > > > mStimulatedAreas
DynamicCellModelLoaderPtr Convert(const FileFinder &rFilePath, bool isCollective=true)
std::string GetMeshName() const
std::vector< std::map< std::string, double > > mParameterSettings
double GetEpiLayerFraction()
virtual double GetParameter(const std::string &rParameterName)=0
double GetDrugDose() const
virtual AbstractLookupTableCollection * GetLookupTableCollection()
AbstractCardiacCellInterface * CreateCardiacCellForTissueNode(Node< SPACE_DIM > *pNode)
std::vector< boost::shared_ptr< AbstractChasteRegion< SPACE_DIM > > > mIonicModelRegions
unsigned GetEpiLayerIndex()
bool IsElectrodesPresent() const
void SetCellIntracellularStimulus(AbstractCardiacCellInterface *pCell, unsigned nodeIndex)
double GetEndoLayerFraction()
unsigned GetMidLayerIndex()
bool GetCheckpointSimulation() const
ChastePoint< SPACE_DIM > GetPoint() const
Definition: Node.cpp:133
unsigned GetIndex() const
Definition: Node.cpp:158
std::vector< cp::ionic_model_selection_type > mIonicModelsDefined
void GetStimuli(std::vector< boost::shared_ptr< AbstractStimulusFunction > > &rStimuliApplied, std::vector< boost::shared_ptr< AbstractChasteRegion< DIM > > > &rStimulatedAreas) const
static HeartConfig * Instance()
void SetCellParameters(AbstractCardiacCellInterface *pCell, unsigned nodeIndex)
std::map< std::string, std::pair< double, double > > GetIc50Values()