BioDynaMo  v1.05.0-137fdb15
cell.h
Go to the documentation of this file.
1 // -----------------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 CERN & University of Surrey for the benefit of the
4 // BioDynaMo collaboration. All Rights Reserved.
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 //
9 // See the LICENSE file distributed with this work for details.
10 // See the NOTICE file distributed with this work for additional information
11 // regarding copyright ownership.
12 //
13 // -----------------------------------------------------------------------------
14 
15 #ifndef CORE_AGENT_CELL_H_
16 #define CORE_AGENT_CELL_H_
17 
18 #include <array>
19 #include <cmath>
20 #include <complex>
21 #include <set>
22 #include <string>
23 #include <type_traits>
24 #include <vector>
25 
26 #include "core/agent/agent.h"
32 #include "core/functor.h"
33 #include "core/interaction_force.h"
34 #include "core/param/param.h"
35 #include "core/shape.h"
36 #include "core/util/math.h"
37 
38 namespace bdm {
39 
40 class Cell : public Agent {
42 
43  public:
45  static const Real3 kXAxis;
47  static const Real3 kYAxis;
49  static const Real3 kZAxis;
50 
51  Cell() : diameter_(1.0), density_(1.0) { UpdateVolume(); }
52 
53  explicit Cell(real_t diameter) : diameter_(diameter), density_(1.0) {
54  UpdateVolume();
55  }
56 
57  explicit Cell(const Real3& position)
58  : position_(position), diameter_(1.0), density_(1.0) {
59  UpdateVolume();
60  }
61 
62  virtual ~Cell() = default;
63 
68  void Initialize(const NewAgentEvent& event) override {
69  Base::Initialize(event);
70 
71  if (event.GetUid() == CellDivisionEvent::kUid) {
72  const auto& cdevent = static_cast<const CellDivisionEvent&>(event);
73  auto* mother_cell = bdm_static_cast<Cell*>(event.existing_agent);
74  auto* daughter = this; // FIXME
75  // A) Defining some values
76  // ..................................................................
77  // defining the two radii s.t total volume is conserved
78  // * radius^3 = r1^3 + r2^3 ;
79  // * volume_ratio = r2^3 / r1^3
80  real_t radius = mother_cell->GetDiameter() * real_t(0.5);
81 
82  // define an axis for division (along which the nuclei will move)
83  real_t x_coord = std::cos(cdevent.theta) * std::sin(cdevent.phi);
84  real_t y_coord = std::sin(cdevent.theta) * std::sin(cdevent.phi);
85  real_t z_coord = std::cos(cdevent.phi);
86  Real3 coords = {x_coord, y_coord, z_coord};
87  real_t total_length_of_displacement = radius / real_t(4.0);
88 
89  const auto& x_axis = mother_cell->kXAxis;
90  const auto& y_axis = mother_cell->kYAxis;
91  const auto& z_axis = mother_cell->kZAxis;
92 
93  Real3 axis_of_division =
94  (coords.EntryWiseProduct(x_axis) + coords.EntryWiseProduct(y_axis) +
95  coords.EntryWiseProduct(z_axis)) *
96  total_length_of_displacement;
97 
98  // two equations for the center displacement :
99  // 1) d2/d1= v2/v1 = volume_ratio (each sphere is shifted inver.
100  // proportionally to its volume)
101  // 2) d1 + d2 = TOTAL_LENGTH_OF_DISPLACEMENT
102  real_t d_2 = total_length_of_displacement / (cdevent.volume_ratio + 1);
103  real_t d_1 = total_length_of_displacement - d_2;
104 
105  real_t mother_volume = mother_cell->GetVolume();
106  real_t new_volume = mother_volume / (cdevent.volume_ratio + 1);
107  daughter->SetVolume(mother_volume - new_volume);
108 
109  // position
110  auto mother_pos = mother_cell->GetPosition();
111  auto new_position = mother_pos + (axis_of_division * d_2);
112  daughter->SetPosition(new_position);
113 
114  // E) This sphere becomes the 1st daughter
115  // move these cells on opposite direction
116  mother_pos -= axis_of_division * d_1;
117  // update mother here and not in Update method to avoid recomputation
118  mother_cell->SetPosition(mother_pos);
119  mother_cell->SetVolume(new_volume);
120 
121  daughter->SetAdherence(mother_cell->GetAdherence());
122  daughter->SetDensity(mother_cell->GetDensity());
123  // G) TODO(lukas) Copy the intracellular and membrane bound Substances
124  }
125  }
126 
127  Shape GetShape() const override { return Shape::kSphere; }
128 
134  virtual Cell* Divide() {
135  auto* random = Simulation::GetActive()->GetRandom();
136  return Divide(random->Uniform(real_t(0.9), real_t(1.1)));
137  }
138 
143  virtual Cell* Divide(real_t volume_ratio) {
144  // find random point on sphere (based on :
145  // http://mathworld.wolfram.com/SpherePointPicking.html)
146  auto* random = Simulation::GetActive()->GetRandom();
147  real_t theta = 2 * Math::kPi * random->Uniform(0, 1);
148  real_t phi = std::acos(2 * random->Uniform(0, 1) - 1);
149  return Divide(volume_ratio, phi, theta);
150  }
151 
156  virtual Cell* Divide(const Real3& axis) {
157  auto* random = Simulation::GetActive()->GetRandom();
158  auto polarcoord = TransformCoordinatesGlobalToPolar(axis + position_);
159  return Divide(random->Uniform(real_t(0.9), real_t(1.1)), polarcoord[1],
160  polarcoord[2]);
161  }
162 
166  virtual Cell* Divide(real_t volume_ratio, const Real3& axis) {
167  auto polarcoord = TransformCoordinatesGlobalToPolar(axis + position_);
168  return Divide(volume_ratio, polarcoord[1], polarcoord[2]);
169  }
170 
174  virtual Cell* Divide(real_t volume_ratio, real_t phi, real_t theta) {
175  CellDivisionEvent event(volume_ratio, phi, theta);
176  CreateNewAgents(event, {this});
177  return bdm_static_cast<Cell*>(event.new_agents[0]);
178  }
179 
180  real_t GetAdherence() const { return adherence_; }
181 
182  real_t GetDiameter() const override { return diameter_; }
183 
184  real_t GetMass() const { return density_ * volume_; }
185 
186  real_t GetDensity() const { return density_; }
187 
188  const Real3& GetPosition() const override { return position_; }
189 
190  const Real3& GetTractorForce() const { return tractor_force_; }
191 
192  real_t GetVolume() const { return volume_; }
193 
194  void SetAdherence(real_t adherence) {
195  if (adherence < adherence_) {
197  }
198  adherence_ = adherence;
199  }
200 
201  void SetDiameter(real_t diameter) override {
202  if (diameter > diameter_) {
204  }
205  diameter_ = diameter;
206  UpdateVolume();
207  }
208 
209  void SetVolume(real_t volume) {
210  volume_ = volume;
211  UpdateDiameter();
212  }
213 
214  void SetMass(real_t mass) { SetDensity(mass / volume_); }
215 
216  void SetDensity(real_t density) {
217  if (density > density_) {
219  }
220  density_ = density;
221  }
222 
223  void SetPosition(const Real3& position) override {
224  position_ = position;
226  }
227 
228  void SetTractorForce(const Real3& tractor_force) {
229  tractor_force_ = tractor_force;
230  }
231 
232  void ChangeVolume(real_t speed) {
233  // scaling for integration step
234  auto* param = Simulation::GetActive()->GetParam();
235  real_t delta = speed * param->simulation_time_step;
236  volume_ += delta;
237  if (volume_ < real_t(5.2359877E-7)) {
238  volume_ = real_t(5.2359877E-7);
239  }
240  UpdateDiameter();
241  }
242 
243  void UpdateDiameter() {
244  // V = (4/3)*pi*r^3 = (pi/6)*diameter^3
245  real_t diameter = std::cbrt(volume_ * 6 / Math::kPi);
246  if (diameter > diameter_) {
247  Base::SetPropagateStaticness();
248  }
249  diameter_ = diameter;
250  }
251 
252  void UpdateVolume() {
253  // V = (4/3)*pi*r^3 = (pi/6)*diameter^3
254  volume_ = Math::kPi / real_t(6) * std::pow(diameter_, 3);
255  }
256 
257  void UpdatePosition(const Real3& delta) {
258  position_ += delta;
260  }
261 
263  real_t squared_radius, real_t dt) override {
264  // Basically, the idea is to make the sum of all the forces acting
265  // on the Point mass. It is stored in translationForceOnPointMass.
266  // There is also a computation of the torque (only applied
267  // by the daughter neurites), stored in rotationForce.
268 
269  // TODO(roman) : There might be a problem, in the sense that the biology
270  // is not applied if the total Force is smaller than adherence.
271  // Once, I should look at this more carefully.
272 
273  // fixme why? copying
274  const auto& tf = GetTractorForce();
275 
276  // the 3 types of movement that can occur
277  // bool biological_translation = false;
278  bool physical_translation = false;
279  // bool physical_rotation = false;
280 
281  real_t h = dt;
282  Real3 movement_at_next_step{0, 0, 0};
283 
284  // BIOLOGY :
285  // 0) Start with tractor force : What the biology defined as active
286  // movement------------
287  movement_at_next_step += tf * h;
288 
289  // PHYSICS
290  // the physics force to move the point mass
291  Real3 translation_force_on_point_mass{0, 0, 0};
292 
293  // the physics force to rotate the cell
294  // Real3 rotation_force { 0, 0, 0 };
295 
296  // 1) "artificial force" to maintain the sphere in the ecm simulation
297  // boundaries--------
298  // 2) Spring force from my neurites (translation and
299  // rotation)--------------------------
300  // 3) Object avoidance force
301  // -----------------------------------------------------------
302  // (We check for every neighbor object if they touch us, i.e. push us
303  // away)
304 
305  uint64_t non_zero_neighbor_forces = 0;
306  if (!IsStatic()) {
307  auto* ctxt = Simulation::GetActive()->GetExecutionContext();
308  auto calculate_neighbor_forces =
309  L2F([&](Agent* neighbor, real_t squared_distance) {
310  auto neighbor_force = force->Calculate(this, neighbor);
311  if (neighbor_force[0] != 0 || neighbor_force[1] != 0 ||
312  neighbor_force[2] != 0) {
313  non_zero_neighbor_forces++;
314  translation_force_on_point_mass[0] += neighbor_force[0];
315  translation_force_on_point_mass[1] += neighbor_force[1];
316  translation_force_on_point_mass[2] += neighbor_force[2];
317  }
318  });
319  ctxt->ForEachNeighbor(calculate_neighbor_forces, *this, squared_radius);
320 
321  if (non_zero_neighbor_forces > 1) {
323  }
324  }
325 
326  // 4) PhysicalBonds
327  // How the physics influences the next displacement
328  real_t norm_of_force = std::sqrt(translation_force_on_point_mass *
329  translation_force_on_point_mass);
330 
331  // is there enough force to :
332  // - make us biologically move (Tractor) :
333  // - break adherence and make us translate ?
334  physical_translation = norm_of_force > GetAdherence();
335 
336  assert(GetMass() != 0 && "The mass of a cell was found to be zero!");
337  real_t mh = h / GetMass();
338  // adding the physics translation (scale by weight) if important enough
339  if (physical_translation) {
340  // We scale the move with mass and time step
341  movement_at_next_step += translation_force_on_point_mass * mh;
342 
343  // Performing the translation itself :
344  // but we want to avoid huge jumps in the simulation, so there are
345  // maximum distances possible
346  auto* param = Simulation::GetActive()->GetParam();
347  if (norm_of_force * mh > param->simulation_max_displacement) {
348  movement_at_next_step.Normalize();
349  movement_at_next_step *= param->simulation_max_displacement;
350  }
351  }
352  return movement_at_next_step;
353  }
354 
355  void ApplyDisplacement(const Real3& displacement) override;
356 
357  void MovePointMass(const Real3& normalized_dir, real_t speed) {
358  tractor_force_ += normalized_dir * speed;
359  }
360 
361  protected:
367  Real3 TransformCoordinatesGlobalToPolar(const Real3& coord) const;
368 
369  private:
371  Real3 position_ = {{0, 0, 0}};
372  Real3 tractor_force_ = {{0, 0, 0}};
380 };
381 
382 } // namespace bdm
383 
384 #endif // CORE_AGENT_CELL_H_
in_place_exec_ctxt.h
bdm::NewAgentEvent
Definition: new_agent_event.h:61
bdm::Agent::SetStaticnessNextTimestep
void SetStaticnessNextTimestep(bool value) const
Definition: agent.h:176
shape.h
bdm::Shape
Shape
Definition: shape.h:20
inline_vector.h
bdm::InteractionForce::Calculate
virtual Real4 Calculate(const Agent *lhs, const Agent *rhs) const
Definition: interaction_force.cc:33
bdm::Cell::density_
real_t density_
NB: Use setter and don't assign values directly.
Definition: cell.h:379
bdm::Cell::Cell
Cell()
Definition: cell.h:51
bdm::Cell::UpdatePosition
void UpdatePosition(const Real3 &delta)
Definition: cell.h:257
bdm
Definition: agent.cc:39
bdm::InteractionForce
Definition: interaction_force.h:26
bdm::Agent::CreateNewAgents
void CreateNewAgents(const NewAgentEvent &event, const std::initializer_list< Agent * > &prototypes)
Definition: agent.h:121
bdm::Cell::Cell
Cell(const Real3 &position)
Definition: cell.h:57
bdm::Cell::BDM_AGENT_HEADER
BDM_AGENT_HEADER(Cell, Agent, 1)
bdm::Cell::GetAdherence
real_t GetAdherence() const
Definition: cell.h:180
bdm::Cell::SetAdherence
void SetAdherence(real_t adherence)
Definition: cell.h:194
bdm::Cell::SetDiameter
void SetDiameter(real_t diameter) override
Definition: cell.h:201
bdm::Cell::Initialize
void Initialize(const NewAgentEvent &event) override
This method is used to initialise the values of daughter 2 for a cell division event.
Definition: cell.h:68
bdm::Cell::diameter_
real_t diameter_
NB: Use setter and don't assign values directly.
Definition: cell.h:374
bdm::Cell::tractor_force_
Real3 tractor_force_
Definition: cell.h:372
bdm::real_t
double real_t
Definition: real_t.h:21
bdm::Cell::UpdateDiameter
void UpdateDiameter()
Definition: cell.h:243
bdm::Cell::Divide
virtual Cell * Divide(real_t volume_ratio, real_t phi, real_t theta)
Divide this cell.
Definition: cell.h:174
bdm::Simulation::GetRandom
Random * GetRandom()
Returns a thread local random number generator.
Definition: simulation.cc:284
bdm::Cell::Cell
Cell(real_t diameter)
Definition: cell.h:53
bdm::Cell::kYAxis
static const Real3 kYAxis
Second axis of the local coordinate system.
Definition: cell.h:47
math.h
bdm::Cell::GetPosition
const Real3 & GetPosition() const override
Definition: cell.h:188
bdm::Agent::SetPropagateStaticness
void SetPropagateStaticness(bool value=true)
Definition: agent.h:184
bdm::Cell::TransformCoordinatesGlobalToPolar
Real3 TransformCoordinatesGlobalToPolar(const Real3 &coord) const
Definition: cell.cc:33
bdm::CellDivisionEvent
Contains the parameters to perform a cell division.
Definition: cell_division_event.h:31
bdm::NewAgentEvent::existing_agent
Agent * existing_agent
Definition: new_agent_event.h:69
bdm::L2F
LambdaFunctor< decltype(&TLambda::operator())> L2F(const TLambda &l)
Definition: functor.h:110
bdm::Agent
Contains code required by all agents.
Definition: agent.h:79
bdm::Cell::GetDensity
real_t GetDensity() const
Definition: cell.h:186
bdm::Cell::Divide
virtual Cell * Divide()
Divide this cell.
Definition: cell.h:134
cell_division_event.h
bdm::Cell::adherence_
real_t adherence_
NB: Use setter and don't assign values directly.
Definition: cell.h:377
bdm::Cell::CalculateDisplacement
Real3 CalculateDisplacement(const InteractionForce *force, real_t squared_radius, real_t dt) override
Definition: cell.h:262
new_agent_event.h
bdm::Cell::GetShape
Shape GetShape() const override
Definition: cell.h:127
bdm::Param::simulation_time_step
real_t simulation_time_step
Definition: param.h:185
bdm::Cell::SetPosition
void SetPosition(const Real3 &position) override
Definition: cell.h:223
bdm::kSphere
@ kSphere
Definition: shape.h:20
math_array.h
bdm::Cell::GetMass
real_t GetMass() const
Definition: cell.h:184
bdm::NewAgentEvent::new_agents
InlineVector< Agent *, 3 > new_agents
Definition: new_agent_event.h:73
bdm::Cell::SetTractorForce
void SetTractorForce(const Real3 &tractor_force)
Definition: cell.h:228
agent.h
bdm::Cell
Definition: cell.h:40
bdm::Cell::ChangeVolume
void ChangeVolume(real_t speed)
Definition: cell.h:232
bdm::Param::simulation_max_displacement
real_t simulation_max_displacement
Definition: param.h:194
bdm::Cell::UpdateVolume
void UpdateVolume()
Definition: cell.h:252
bdm::Cell::GetDiameter
real_t GetDiameter() const override
Definition: cell.h:182
bdm::Cell::SetDensity
void SetDensity(real_t density)
Definition: cell.h:216
bdm::Simulation::GetExecutionContext
ExecutionContext * GetExecutionContext()
Returns a thread local execution context.
Definition: simulation.cc:288
bdm::Cell::GetTractorForce
const Real3 & GetTractorForce() const
Definition: cell.h:190
bdm::Cell::SetMass
void SetMass(real_t mass)
Definition: cell.h:214
bdm::Cell::GetVolume
real_t GetVolume() const
Definition: cell.h:192
bdm::NewAgentEvent::GetUid
virtual NewAgentEventUid GetUid() const =0
bdm::Cell::SetVolume
void SetVolume(real_t volume)
Definition: cell.h:209
bdm::Cell::kXAxis
static const Real3 kXAxis
First axis of the local coordinate system.
Definition: cell.h:45
bdm::Math::kPi
static constexpr real_t kPi
value of pi
Definition: math.h:33
bdm::MathArray::EntryWiseProduct
MathArray EntryWiseProduct(const MathArray &rhs) const
Definition: math_array.h:396
bdm::Simulation::GetParam
const Param * GetParam() const
Returns the simulation parameters.
Definition: simulation.cc:271
bdm::Cell::Divide
virtual Cell * Divide(const Real3 &axis)
Divide this cell.
Definition: cell.h:156
bdm::MathArray< real_t, 3 >
bdm::Cell::MovePointMass
void MovePointMass(const Real3 &normalized_dir, real_t speed)
Definition: cell.h:357
bdm::Cell::ApplyDisplacement
void ApplyDisplacement(const Real3 &displacement) override
Definition: cell.cc:23
bdm::Cell::~Cell
virtual ~Cell()=default
param.h
bdm::Simulation::GetActive
static Simulation * GetActive()
This function returns the currently active Simulation simulation.
Definition: simulation.cc:67
bdm::Cell::Divide
virtual Cell * Divide(real_t volume_ratio)
Divide this cell.
Definition: cell.h:143
bdm::Cell::Divide
virtual Cell * Divide(real_t volume_ratio, const Real3 &axis)
Divide this cell.
Definition: cell.h:166
bdm::CellDivisionEvent::kUid
static const NewAgentEventUid kUid
Definition: cell_division_event.h:32
functor.h
bdm::Agent::IsStatic
bool IsStatic() const
Definition: agent.h:204
interaction_force.h
bdm::Cell::kZAxis
static const Real3 kZAxis
Third axis of the local coordinate system.
Definition: cell.h:49
bdm::Cell::position_
Real3 position_
NB: Use setter and don't assign values directly.
Definition: cell.h:371
bdm::Cell::volume_
real_t volume_
Definition: cell.h:375