Demo monolayer growth

demo/monolayer_growth/CMakeLists.txt

@@ -0,0 +1,35 @@
+# -----------------------------------------------------------------------------
+#
+# Copyright (C) 2022 CERN & University of Surrey for the benefit of the
+# BioDynaMo collaboration. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+#
+# See the LICENSE file distributed with this work for details.
+# See the NOTICE file distributed with this work for additional information
+# regarding copyright ownership.
+#
+# -----------------------------------------------------------------------------
+cmake_minimum_required(VERSION 3.19.3)
+project(monolayer_growth)
+
+# BioDynaMo curretly uses the C++14 standard.
+set(CMAKE_CXX_STANDARD 14)
+
+# Use BioDynaMo in this project.
+find_package(BioDynaMo REQUIRED)
+
+# See UseBioDynaMo.cmake in your BioDynaMo build folder for details.
+# Note that BioDynaMo provides gtest header/libraries in its include/lib dir.
+include(${BDM_USE_FILE})
+
+# Consider all files in src/ for BioDynaMo simulation.
+include_directories("src")
+file(GLOB_RECURSE PROJECT_HEADERS src/*.h)
+file(GLOB_RECURSE PROJECT_SOURCES src/*.cc)
+
+bdm_add_executable(${CMAKE_PROJECT_NAME}
+                   HEADERS ${PROJECT_HEADERS}
+                   SOURCES ${PROJECT_SOURCES}
+                   LIBRARIES ${BDM_REQUIRED_LIBRARIES})

demo/monolayer_growth/src/behaviours.h

@@ -0,0 +1,116 @@
+// -----------------------------------------------------------------------------
+//
+// Copyright (C) 2021 CERN & University of Surrey for the benefit of the
+// BioDynaMo collaboration. All Rights Reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+//
+// See the LICENSE file distributed with this work for details.
+// See the NOTICE file distributed with this work for additional information
+// regarding copyright ownership.
+//
+// -----------------------------------------------------------------------------
+
+#ifndef BEHAVIOURS_H_
+#define BEHAVIOURS_H_
+
+#include "biodynamo.h"
+#include "core/behavior/behavior.h"
+#include "cycling_cell.h"
+#include "sim_param.h"
+
+namespace bdm {
+
+// Define growth behaviour
+struct GrowthAndCellCycle : public Behavior {
+  BDM_BEHAVIOR_HEADER(GrowthAndCellCycle, Behavior, 1);
+
+  GrowthAndCellCycle() { AlwaysCopyToNew(); }
+  virtual ~GrowthAndCellCycle() {}
+
+  void Run(Agent* agent) override {
+    // Define the maximum times for each phase of the cell cycle
+    const real_t kG1Duration = 7;
+    const real_t kSDuration = 6;
+    const real_t kG2Duration = 3;
+    const real_t kMDuration = 2;
+
+    if (auto* cell = dynamic_cast<CyclingCell*>(agent)) {
+      auto* random = Simulation::GetActive()->GetRandom();
+      real_t ran = random->Uniform(0, 1) * 1.0;
+      const auto* param = Simulation::GetActive()->GetParam();
+      const auto* sparam = param->Get<SimParam>();
+
+      // Counter for Delta t at each stage
+      // Used for calculating probability of moving to next state.
+      cell->SetDeltaT(cell->GetDeltaT() + param->simulation_time_step);
+
+      // If statements for checking what states we are in and if
+      // a cell moves to the next state based on cumulative probability.
+      if (cell->GetCycle() == CellState::kG1) {
+        real_t p1 = (cell->GetDeltaT() / kG1Duration);
+        if (p1 > ran) {
+          // Changing cells state number or "cycle" postiion.
+          cell->SetCycle(CellState::kS);
+          // Delta t is always reset when exiting a state for use in the next
+          // state.
+          cell->SetDeltaT(0);
+        }
+
+      } else if (cell->GetCycle() == CellState::kS) {
+        real_t p2 = (cell->GetDeltaT() / kSDuration);
+        if (p2 > ran) {
+          cell->SetCycle(CellState::kG2);
+          cell->SetDeltaT(0);
+        }
+
+      } else if (cell->GetCycle() == CellState::kG2) {
+        real_t p3 = (cell->GetDeltaT() / kG2Duration);
+        if (p3 > ran) {
+          cell->SetCycle(CellState::kM);
+          cell->SetDeltaT(0);
+        }
+
+      } else {
+        real_t p4 = (cell->GetDeltaT() / kMDuration);
+        if (p4 > ran) {
+          cell->SetCycle(CellState::kG1);
+          cell->SetDeltaT(0);
+          // Checking if cell has reached the critical volume which leads to
+          // cell division. Here 0.975 Vmax is roughly 195% the initial cell
+          // volume.
+          if (cell->GetVolume() > cell->GetVmax() * 0.975) {
+            int neighbours_counter = 0;
+            auto* ctxt = Simulation::GetActive()->GetExecutionContext();
+            auto check_surrounding_cells =
+                L2F([&](Agent* neighbor, real_t squared_distance) {
+                  neighbours_counter++;
+                });
+            ctxt->ForEachNeighbor(check_surrounding_cells, *cell,
+                                  (pow(cell->GetDiameter(), 2)));
+            // Divide only if the number  of neighbours doesn't exceed a
+            // predefined threshold (might represent direct contact signal,
+            // pressure, etc)
+            if (neighbours_counter <= sparam->neighbours_threshold) {
+              cell->Divide(1);
+            }
+          }
+        }
+      }
+
+      // Checking if our cells volume is less than the maximum possible valuable
+      // achievalbe. if yes cell grows if no then cell does not grow.
+      if (cell->GetVolume() < cell->GetVmax()) {
+        real_t alpha = 1.0;
+        cell->ChangeVolume(
+            alpha * cell->GetVolume() *
+            ((cell->GetVmax() - cell->GetVolume()) / cell->GetVmax()));
+      }
+    }
+  }
+};
+
+}  // namespace bdm
+
+#endif  // BEHAVIOURS_H_

demo/monolayer_growth/src/cell_cell_force.cc

@@ -0,0 +1,75 @@
+// -----------------------------------------------------------------------------
+//
+// Copyright (C) 2021 CERN & University of Surrey for the benefit of the
+// BioDynaMo collaboration. All Rights Reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+//
+// See the LICENSE file distributed with this work for details.
+// See the NOTICE file distributed with this work for additional information
+// regarding copyright ownership.
+//
+// -----------------------------------------------------------------------------
+
+#include "cell_cell_force.h"
+#include "sim_param.h"
+
+namespace bdm {
+
+/// Custom force. Changed adhesive and repulsive parameters compared to standard
+/// force to achieve quick separation of mother and daughter cells after
+/// division.
+Real4 CellCellForce::Calculate(const Agent* lhs, const Agent* rhs) const {
+  const Real3& ref_mass_location = lhs->GetPosition();
+  real_t ref_diameter = lhs->GetDiameter();
+  real_t ref_iof_coefficient = 0.15;
+  const Real3& nb_mass_location = rhs->GetPosition();
+  real_t nb_diameter = rhs->GetDiameter();
+  real_t nb_iof_coefficient = 0.15;
+
+  auto c1 = ref_mass_location;
+  real_t r1 = 0.5 * ref_diameter;
+  auto c2 = nb_mass_location;
+  real_t r2 = 0.5 * nb_diameter;
+  // We take virtual bigger radii to have a distant interaction, to get a
+  // desired density.
+  real_t additional_radius =
+      10.0 * std::min(ref_iof_coefficient, nb_iof_coefficient);
+  r1 += additional_radius;
+  r2 += additional_radius;
+  // the 3 components of the vector c2 -> c1
+  real_t comp1 = c1[0] - c2[0];
+  real_t comp2 = c1[1] - c2[1];
+  real_t comp3 = c1[2] - c2[2];
+  real_t center_distance =
+      std::sqrt(comp1 * comp1 + comp2 * comp2 + comp3 * comp3);
+  // the overlap distance (how much one penetrates in the other)
+  real_t delta = r1 + r2 - center_distance;
+  // if no overlap : no force
+  if (delta < 0) {
+    return {0.0, 0.0, 0.0, 0.0};
+  }
+  // to avoid a division by 0 if the centers are (almost) at the same
+  //  location
+  if (center_distance < 0.00000001) {
+    auto* random = Simulation::GetActive()->GetRandom();
+    auto force2on1 = random->template UniformArray<3>(-3.0, 3.0);
+    return {force2on1[0], force2on1[1], force2on1[2], 0};
+  }
+  // the force itself
+  const auto* sparam =
+      Simulation::GetActive()
+          ->GetParam()
+          ->Get<SimParam>();  // get a pointer to an instance of SimParam
+  real_t r = (r1 * r2) / (r1 + r2);
+  real_t gamma = sparam->attraction_coeff;
+  real_t k = sparam->repulsion_coeff;
+  real_t f = k * delta - gamma * std::sqrt(r * delta);
+
+  real_t module = f / center_distance;
+  Real3 force2on1({module * comp1, module * comp2, module * comp3});
+  return {force2on1[0], force2on1[1], force2on1[2], 0};
+}
+
+}  // namespace bdm

demo/monolayer_growth/src/cell_cell_force.h

@@ -0,0 +1,33 @@
+// -----------------------------------------------------------------------------
+//
+// Copyright (C) 2021 CERN & University of Surrey for the benefit of the
+// BioDynaMo collaboration. All Rights Reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+//
+// See the LICENSE file distributed with this work for details.
+// See the NOTICE file distributed with this work for additional information
+// regarding copyright ownership.
+//
+// -----------------------------------------------------------------------------
+
+#ifndef CELL_CELL_FORCE_H_
+#define CELL_CELL_FORCE_H_
+
+#include "biodynamo.h"
+#include "core/interaction_force.h"
+
+namespace bdm {
+
+class CellCellForce : public InteractionForce {
+ public:
+  CellCellForce() {}
+  virtual ~CellCellForce() {}
+
+  virtual Real4 Calculate(const Agent* lhs, const Agent* rhs) const override;
+};
+
+}  // namespace bdm
+
+#endif

demo/monolayer_growth/src/custom_ops.cc

@@ -0,0 +1,19 @@
+// -----------------------------------------------------------------------------
+//
+// Copyright (C) 2021 CERN & University of Surrey for the benefit of the
+// BioDynaMo collaboration. All Rights Reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+//
+// See the LICENSE file distributed with this work for details.
+// See the NOTICE file distributed with this work for additional information
+// regarding copyright ownership.
+//
+// -----------------------------------------------------------------------------
+
+#include "custom_ops.h"
+
+namespace bdm {
+BDM_REGISTER_OP(MoveCellsPlane, "move_cells_plane", kCpu);
+}  // namespace bdm

demo/monolayer_growth/src/custom_ops.h

@@ -0,0 +1,38 @@
+// -----------------------------------------------------------------------------
+//
+// Copyright (C) 2021 CERN & University of Surrey for the benefit of the
+// BioDynaMo collaboration. All Rights Reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+//
+// See the LICENSE file distributed with this work for details.
+// See the NOTICE file distributed with this work for additional information
+// regarding copyright ownership.
+//
+// -----------------------------------------------------------------------------
+
+#ifndef CUSTOM_OPS_H_
+#define CUSTOM_OPS_H_
+
+#include "biodynamo.h"
+#include "core/environment/uniform_grid_environment.h"
+
+namespace bdm {
+
+// Brings cells back to the xy plane (i.e. z = 0)
+struct MoveCellsPlane : public AgentOperationImpl {
+  BDM_OP_HEADER(MoveCellsPlane);
+
+  void operator()(Agent* agent) override {
+    Real3 cell_position = agent->GetPosition();
+    if (cell_position[2] != 0.0) {
+      cell_position[2] = 0.0;
+      agent->SetPosition(cell_position);
+    }
+  }
+};
+
+}  // namespace bdm
+
+#endif