Chaste  Release::2018.1
NagaiHondaForce.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 "NagaiHondaForce.hpp"
37 
38 template<unsigned DIM>
40  : AbstractForce<DIM>(),
41  mNagaiHondaDeformationEnergyParameter(100.0), // This is 1.0 in the Nagai & Honda paper.
42  mNagaiHondaMembraneSurfaceEnergyParameter(10.0), // This is 0.1 in the Nagai & Honda paper.
43  mNagaiHondaCellCellAdhesionEnergyParameter(0.5), // This corresponds to a value of 1.0 for
44  // the sigma parameter in the Nagai & Honda
45  // paper. In the paper, the sigma value is
46  // set to 0.01.
47  mNagaiHondaCellBoundaryAdhesionEnergyParameter(1.0) // This is 0.01 in the Nagai & Honda paper.
48 {
49 }
50 
51 template<unsigned DIM>
53 {
54 }
55 
56 template<unsigned DIM>
58 {
59  // Throw an exception message if not using a VertexBasedCellPopulation
60  if (dynamic_cast<VertexBasedCellPopulation<DIM>*>(&rCellPopulation) == nullptr)
61  {
62  EXCEPTION("NagaiHondaForce is to be used with a VertexBasedCellPopulation only");
63  }
64 
65  // Define some helper variables
66  VertexBasedCellPopulation<DIM>* p_cell_population = static_cast<VertexBasedCellPopulation<DIM>*>(&rCellPopulation);
67  unsigned num_nodes = p_cell_population->GetNumNodes();
68  unsigned num_elements = p_cell_population->GetNumElements();
69 
70  // Begin by computing the area and perimeter of each element in the mesh, to avoid having to do this multiple times
71  std::vector<double> element_areas(num_elements);
72  std::vector<double> element_perimeters(num_elements);
73  std::vector<double> target_areas(num_elements);
74  for (typename VertexMesh<DIM,DIM>::VertexElementIterator elem_iter = p_cell_population->rGetMesh().GetElementIteratorBegin();
75  elem_iter != p_cell_population->rGetMesh().GetElementIteratorEnd();
76  ++elem_iter)
77  {
78  unsigned elem_index = elem_iter->GetIndex();
79  element_areas[elem_index] = p_cell_population->rGetMesh().GetVolumeOfElement(elem_index);
80  element_perimeters[elem_index] = p_cell_population->rGetMesh().GetSurfaceAreaOfElement(elem_index);
81  try
82  {
83  // If we haven't specified a growth modifier, there won't be any target areas in the CellData array and CellData
84  // will throw an exception that it doesn't have "target area" entries. We add this piece of code to give a more
85  // understandable message. There is a slight chance that the exception is thrown although the error is not about the
86  // target areas.
87  target_areas[elem_index] = p_cell_population->GetCellUsingLocationIndex(elem_index)->GetCellData()->GetItem("target area");
88  }
89  catch (Exception&)
90  {
91  EXCEPTION("You need to add an AbstractTargetAreaModifier to the simulation in order to use NagaiHondaForce");
92  }
93  }
94 
95  // Iterate over vertices in the cell population
96  for (unsigned node_index=0; node_index<num_nodes; node_index++)
97  {
98  Node<DIM>* p_this_node = p_cell_population->GetNode(node_index);
99 
100  /*
101  * The force on this Node is given by the gradient of the total free
102  * energy of the CellPopulation, evaluated at the position of the vertex. This
103  * free energy is the sum of the free energies of all CellPtrs in
104  * the cell population. The free energy of each CellPtr is comprised of three
105  * parts - a cell deformation energy, a membrane surface tension energy
106  * and an adhesion energy.
107  *
108  * Note that since the movement of this Node only affects the free energy
109  * of the CellPtrs containing it, we can just consider the contributions
110  * to the free energy gradient from each of these CellPtrs.
111  */
112  c_vector<double, DIM> deformation_contribution = zero_vector<double>(DIM);
113  c_vector<double, DIM> membrane_surface_tension_contribution = zero_vector<double>(DIM);
114  c_vector<double, DIM> adhesion_contribution = zero_vector<double>(DIM);
115 
116  // Find the indices of the elements owned by this node
117  std::set<unsigned> containing_elem_indices = p_cell_population->GetNode(node_index)->rGetContainingElementIndices();
118 
119  // Iterate over these elements
120  for (std::set<unsigned>::iterator iter = containing_elem_indices.begin();
121  iter != containing_elem_indices.end();
122  ++iter)
123  {
124  // Get this element, its index and its number of nodes
125  VertexElement<DIM, DIM>* p_element = p_cell_population->GetElement(*iter);
126  unsigned elem_index = p_element->GetIndex();
127  unsigned num_nodes_elem = p_element->GetNumNodes();
128 
129  // Find the local index of this node in this element
130  unsigned local_index = p_element->GetNodeLocalIndex(node_index);
131 
132  // Add the force contribution from this cell's deformation energy (note the minus sign)
133  c_vector<double, DIM> element_area_gradient = p_cell_population->rGetMesh().GetAreaGradientOfElementAtNode(p_element, local_index);
134  deformation_contribution -= 2*GetNagaiHondaDeformationEnergyParameter()*(element_areas[elem_index] - target_areas[elem_index])*element_area_gradient;
135 
136  // Get the previous and next nodes in this element
137  unsigned previous_node_local_index = (num_nodes_elem+local_index-1)%num_nodes_elem;
138  Node<DIM>* p_previous_node = p_element->GetNode(previous_node_local_index);
139 
140  unsigned next_node_local_index = (local_index+1)%num_nodes_elem;
141  Node<DIM>* p_next_node = p_element->GetNode(next_node_local_index);
142 
143  // Compute the adhesion parameter for each of these edges
144  double previous_edge_adhesion_parameter = GetAdhesionParameter(p_previous_node, p_this_node, *p_cell_population);
145  double next_edge_adhesion_parameter = GetAdhesionParameter(p_this_node, p_next_node, *p_cell_population);
146 
147  // Compute the gradient of each these edges, computed at the present node
148  c_vector<double, DIM> previous_edge_gradient = -p_cell_population->rGetMesh().GetNextEdgeGradientOfElementAtNode(p_element, previous_node_local_index);
149  c_vector<double, DIM> next_edge_gradient = p_cell_population->rGetMesh().GetNextEdgeGradientOfElementAtNode(p_element, local_index);
150 
151  // Add the force contribution from cell-cell and cell-boundary adhesion (note the minus sign)
152  adhesion_contribution -= previous_edge_adhesion_parameter*previous_edge_gradient + next_edge_adhesion_parameter*next_edge_gradient;
153 
154  // Add the force contribution from this cell's membrane surface tension (note the minus sign)
155  c_vector<double, DIM> element_perimeter_gradient = previous_edge_gradient + next_edge_gradient;
156  double cell_target_perimeter = 2*sqrt(M_PI*target_areas[elem_index]);
157  membrane_surface_tension_contribution -= 2*GetNagaiHondaMembraneSurfaceEnergyParameter()*(element_perimeters[elem_index] - cell_target_perimeter)*element_perimeter_gradient;
158  }
159 
160  c_vector<double, DIM> force_on_node = deformation_contribution + membrane_surface_tension_contribution + adhesion_contribution;
161  p_cell_population->GetNode(node_index)->AddAppliedForceContribution(force_on_node);
162  }
163 }
164 
165 template<unsigned DIM>
167 {
168  // Find the indices of the elements owned by each node
169  std::set<unsigned> elements_containing_nodeA = pNodeA->rGetContainingElementIndices();
170  std::set<unsigned> elements_containing_nodeB = pNodeB->rGetContainingElementIndices();
171 
172  // Find common elements
173  std::set<unsigned> shared_elements;
174  std::set_intersection(elements_containing_nodeA.begin(),
175  elements_containing_nodeA.end(),
176  elements_containing_nodeB.begin(),
177  elements_containing_nodeB.end(),
178  std::inserter(shared_elements, shared_elements.begin()));
179 
180  // Check that the nodes have a common edge
181  assert(!shared_elements.empty());
182 
183  double adhesion_parameter = GetNagaiHondaCellCellAdhesionEnergyParameter();
184 
185  // If the edge corresponds to a single element, then the cell is on the boundary
186  if (shared_elements.size() == 1)
187  {
188  adhesion_parameter = GetNagaiHondaCellBoundaryAdhesionEnergyParameter();
189  }
190 
191  return adhesion_parameter;
192 }
193 
194 template<unsigned DIM>
196 {
197  return mNagaiHondaDeformationEnergyParameter;
198 }
199 
200 template<unsigned DIM>
202 {
203  return mNagaiHondaMembraneSurfaceEnergyParameter;
204 }
205 
206 template<unsigned DIM>
208 {
209  return mNagaiHondaCellCellAdhesionEnergyParameter;
210 }
211 
212 template<unsigned DIM>
214 {
215  return mNagaiHondaCellBoundaryAdhesionEnergyParameter;
216 }
217 
218 template<unsigned DIM>
220 {
221  mNagaiHondaDeformationEnergyParameter = deformationEnergyParameter;
222 }
223 
224 template<unsigned DIM>
225 void NagaiHondaForce<DIM>::SetNagaiHondaMembraneSurfaceEnergyParameter(double membraneSurfaceEnergyParameter)
226 {
227  mNagaiHondaMembraneSurfaceEnergyParameter = membraneSurfaceEnergyParameter;
228 }
229 
230 template<unsigned DIM>
231 void NagaiHondaForce<DIM>::SetNagaiHondaCellCellAdhesionEnergyParameter(double cellCellAdhesionEnergyParameter)
232 {
233  mNagaiHondaCellCellAdhesionEnergyParameter = cellCellAdhesionEnergyParameter;
234 }
235 
236 template<unsigned DIM>
237 void NagaiHondaForce<DIM>::SetNagaiHondaCellBoundaryAdhesionEnergyParameter(double cellBoundaryAdhesionEnergyParameter)
238 {
239  mNagaiHondaCellBoundaryAdhesionEnergyParameter = cellBoundaryAdhesionEnergyParameter;
240 }
241 
242 template<unsigned DIM>
243 void NagaiHondaForce<DIM>::OutputForceParameters(out_stream& rParamsFile)
244 {
245  *rParamsFile << "\t\t\t<NagaiHondaDeformationEnergyParameter>" << mNagaiHondaDeformationEnergyParameter << "</NagaiHondaDeformationEnergyParameter>\n";
246  *rParamsFile << "\t\t\t<NagaiHondaMembraneSurfaceEnergyParameter>" << mNagaiHondaMembraneSurfaceEnergyParameter << "</NagaiHondaMembraneSurfaceEnergyParameter>\n";
247  *rParamsFile << "\t\t\t<NagaiHondaCellCellAdhesionEnergyParameter>" << mNagaiHondaCellCellAdhesionEnergyParameter << "</NagaiHondaCellCellAdhesionEnergyParameter>\n";
248  *rParamsFile << "\t\t\t<NagaiHondaCellBoundaryAdhesionEnergyParameter>" << mNagaiHondaCellBoundaryAdhesionEnergyParameter << "</NagaiHondaCellBoundaryAdhesionEnergyParameter>\n";
249 
250  // Call method on direct parent class
252 }
253 
254 // Explicit instantiation
255 template class NagaiHondaForce<1>;
256 template class NagaiHondaForce<2>;
257 template class NagaiHondaForce<3>;
258 
259 // Serialization for Boost >= 1.36
virtual CellPtr GetCellUsingLocationIndex(unsigned index)
Definition: Node.hpp:58
#define EXCEPTION(message)
Definition: Exception.hpp:143
std::set< unsigned > & rGetContainingElementIndices()
Definition: Node.cpp:300
c_vector< double, SPACE_DIM > GetNextEdgeGradientOfElementAtNode(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned localIndex)
double GetNagaiHondaMembraneSurfaceEnergyParameter()
void SetNagaiHondaMembraneSurfaceEnergyParameter(double nagaiHondaMembraneSurfaceEnergyParameter)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
MutableVertexMesh< DIM, DIM > & rGetMesh()
VertexElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
Definition: VertexMesh.hpp:679
void SetNagaiHondaCellBoundaryAdhesionEnergyParameter(double nagaiHondaCellBoundaryAdhesionEnergyParameter)
void OutputForceParameters(out_stream &rParamsFile)
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
Node< DIM > * GetNode(unsigned index)
unsigned GetNumNodes() const
void SetNagaiHondaDeformationEnergyParameter(double nagaiHondaDeformationEnergyParameter)
virtual ~NagaiHondaForce()
virtual double GetVolumeOfElement(unsigned index)
void SetNagaiHondaCellCellAdhesionEnergyParameter(double nagaiHondaCellCellAdhesionEnergyEnergyParameter)
double GetNagaiHondaCellCellAdhesionEnergyParameter()
unsigned GetIndex() const
double GetNagaiHondaCellBoundaryAdhesionEnergyParameter()
double GetNagaiHondaDeformationEnergyParameter()
c_vector< double, SPACE_DIM > GetAreaGradientOfElementAtNode(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned localIndex)
unsigned GetNodeLocalIndex(unsigned globalIndex) const
virtual double GetAdhesionParameter(Node< DIM > *pNodeA, Node< DIM > *pNodeB, VertexBasedCellPopulation< DIM > &rVertexCellPopulation)
VertexElementIterator GetElementIteratorEnd()
Definition: VertexMesh.hpp:686
virtual double GetSurfaceAreaOfElement(unsigned index)
VertexElement< DIM, DIM > * GetElement(unsigned elementIndex)
virtual void AddForceContribution(AbstractCellPopulation< DIM > &rCellPopulation)
virtual void OutputForceParameters(out_stream &rParamsFile)=0