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