Model with CAR-T cells

.gitignore

@@ -121,4 +121,7 @@ doc/latex
 *.dat
 
 # Draft code for comparing models
-draft_code_my_own_analysis/
+draft_code_my_own_analysis/
+
+# Temporary files
+cout.txt

CMakeLists.txt

@@ -1,5 +1,5 @@
 cmake_minimum_required(VERSION 3.19.3)
-project(cart_tumor)
+project(CARTopiaX)
 
 
 # BioDynaMo curretly uses the C++17 standard.

bdm.toml

@@ -10,5 +10,12 @@ additional_data_members = [
 	{name = "cell_type", function = "GetTypeAsInt"}
 ]
 
+[[visualize_agent]]
+name = "CarTCell"
+additional_data_members = ["diameter_","volume_"]
+
 [[visualize_diffusion]]
 name = "oxygen"
+
+[[visualize_diffusion]]
+name = "immunostimulatory_factor"

src/cart_cell.cc

@@ -24,14 +24,17 @@
 #include "tumor_cell.h"
 #include "utils_aux.h"
 #include "core/agent/agent.h"
+#include "core/agent/agent_pointer.h"
 #include "core/agent/new_agent_event.h"
 #include "core/container/math_array.h"
 #include "core/diffusion/diffusion_grid.h"
 #include "core/functor.h"
 #include "core/interaction_force.h"
 #include "core/real_t.h"
 #include "core/resource_manager.h"
+#include "core/simulation.h"
 #include "core/util/log.h"
+#include "core/util/random.h"
 #include <algorithm>
 #include <cmath>
 #include <cstdint>
@@ -47,6 +50,12 @@ CarTCell::CarTCell(const Real3& position) {
   oxygen_dgrid_ = rm.GetDiffusionGrid("oxygen");
   immunostimulatory_factor_dgrid_ =
       rm.GetDiffusionGrid("immunostimulatory_factor");
+
+  SetCurrentLiveTime(kAverageMaximumTimeUntillApoptosisCart);
+  // Add Consumption and Secretion
+  //  Set default oxygen consumption rate
+  SetOxygenConsumptionRate(kDefaultOxygenConsumptionCarT);
+  // Compute constants for all ConsumptionSecretion of Oxygen
   ComputeConstantsConsumptionSecretion();
 }
 
@@ -138,6 +147,8 @@ void CarTCell::ChangeVolumeExponentialRelaxationEquation(
 // compute Displacement
 Real3 CarTCell::CalculateDisplacement(const InteractionForce* force,
                                       real_t squared_radius, real_t /*dt*/) {
+  Simulation* sim = Simulation::GetActive();
+
   // real_t h = dt;
   Real3 movement_at_next_step{0, 0, 0};
   // this should be chaged in a future version of BioDynaMo in order to have a
@@ -148,33 +159,114 @@ Real3 CarTCell::CalculateDisplacement(const InteractionForce* force,
 
   Real3 translation_velocity_on_point_mass{0, 0, 0};
 
-  // We check for every neighbor object if they touch us, i.e. push us
-  // away and agreagate the velocities
+  //--------------------------------------------
+  // CAR-T self motility (in case of migration)
+  //--------------------------------------------
+  Real3 current_position = GetPosition();
+  ExecutionContext* ctxt = sim->GetExecutionContext();
+  Random* rng = sim->GetRandom();
+  Real3 motility;
+  if (DoesCellMove()) {
+    // compute motility
+    if (rng->Uniform(0.0, 1.0) < kMotilityProbability) {
+      // random direction as unitary vector
+      const Real3 random_direction = GenerateRandomDirection();
+      Real3 direction_to_immunostimulatory_factor;
+      // returns normalized gradient towards the immunostimulatory factor source
+      immunostimulatory_factor_dgrid_->GetGradient(
+          current_position, &direction_to_immunostimulatory_factor, true);
+      // motility = bias * direction_to_immunostimulatory_factor +
+      // (1-bias)*random_direction
+      motility = kMigrationBiasCart * direction_to_immunostimulatory_factor +
+                 kMigrationOneMinusBiasCart * random_direction;
+      const real_t motility_norm_squared = motility[0] * motility[0] +
+                                           motility[1] * motility[1] +
+                                           motility[2] * motility[2];
+      // Convert to unit direction
+      if (motility_norm_squared > 0) {
+        motility.Normalize();
+      }
+      // Scale by migration speed and add to the velocity
+      translation_velocity_on_point_mass += motility * kMigrationSpeedCart;
+    }
+  }
+
+  //--------------------------------------------
+  // CAR-T killing or victim cell escaping
+  //--------------------------------------------
+  // If cell is not apoptotic
+  if (state_ == CarTCellState::kAlive) {
+    // if it is attached to tumor cell
+    if (attached_to_tumor_cell_) {
+      // try to kill the cancer cell and in case of failure see if it manages to
+      // scape the order needs to be this one because it should try to kill
+      // before seeing if it scapes
+      if (TryToInduceApoptosis(attached_cell_ptr_, rng) ||
+          rng->Uniform(0.0, 1.0) < kProbabilityEscape) {
+        // the cancer cell is detached
+        attached_cell_ptr_->SetAttachedToCart(false);
+        // empty ID
+        attached_cell_ptr_ = nullptr;
+        attached_to_tumor_cell_ = false;
+      }
+    }
 
-  uint64_t non_zero_neighbor_forces = 0;
-  if (!IsStatic()) {
-    auto* ctxt = Simulation::GetActive()->GetExecutionContext();
-    auto calculate_neighbor_forces =
+    //--------------------------------------------
+    // CAR-T adhesion to victim cell
+    //--------------------------------------------
+    // Compute forces between the cells and check for a new attachment
+    auto calculate_forces_and_elastic_displacement =
         L2F([&](Agent* neighbor, real_t /*squared_distance*/) {
-          auto neighbor_force = force->Calculate(this, neighbor);
-          if (neighbor_force[0] != 0 || neighbor_force[1] != 0 ||
-              neighbor_force[2] != 0) {
-            non_zero_neighbor_forces++;
-            translation_velocity_on_point_mass[0] += neighbor_force[0];
-            translation_velocity_on_point_mass[1] += neighbor_force[1];
-            translation_velocity_on_point_mass[2] += neighbor_force[2];
+          // Adhesion repulsion forces between cells
+          //  We check for every neighbor object if they touch us, i.e. push us
+          //  away and aggregate the velocities
+          Real4 neighbor_force = force->Calculate(this, neighbor);
+          translation_velocity_on_point_mass[0] += neighbor_force[0];
+          translation_velocity_on_point_mass[1] += neighbor_force[1];
+          translation_velocity_on_point_mass[2] += neighbor_force[2];
+
+          // CAR-T adhesion to new victim cell
+          Real3 displac = neighbor->GetPosition() - current_position;
+
+          if (auto* cancer_cell = dynamic_cast<TumorCell*>(neighbor)) {
+            // movement towards the tumor cells
+            const real_t sq_norm_displac = displac[0] * displac[0] +
+                                           displac[1] * displac[1] +
+                                           displac[2] * displac[2];
+
+            // The cart moves towards the tumor cell only if they are not
+            // touching already If they are too close the only force affecting
+            // is the adhesion force to avoid CAR-T non-stop pushing tumor
+            // cells. In case of being closer than
+            // gKMaxSquaredDistanceCartMovingTowardsTumorCell there is a
+            // probability kProbabilityPushing for the CAR-T to keep pushing the
+            // tumor cell
+            if (sq_norm_displac >
+                gKMaxSquaredDistanceCartMovingTowardsTumorCell) {
+              translation_velocity_on_point_mass[0] +=
+                  displac[0] * kElasticConstantCart;
+              translation_velocity_on_point_mass[1] +=
+                  displac[1] * kElasticConstantCart;
+              translation_velocity_on_point_mass[2] +=
+                  displac[2] * kElasticConstantCart;
+            }
+
+            // If the CAR-T has not succeeded in attaching to a tumor cell yet,
+            // it tries again
+            if (!attached_to_tumor_cell_) {
+              TryToGetAttachedTo(cancer_cell, sq_norm_displac, rng);
+            }
           }
         });
-    ctxt->ForEachNeighbor(calculate_neighbor_forces, *this, squared_radius);
-
-    if (non_zero_neighbor_forces > 1) {
-      SetStaticnessNextTimestep(false);
-    }
+    ctxt->ForEachNeighbor(calculate_forces_and_elastic_displacement, *this,
+                          squared_radius);
   }
 
+  //--------------------------------------------
   // Two step Adams-Bashforth approximation of the time derivative for position
   // position(t + dt) ≈ position(t) + dt * [ 1.5 * velocity(t) - 0.5 *
   // velocity(t - dt) ]
+  //--------------------------------------------
   movement_at_next_step +=
       translation_velocity_on_point_mass * kDnew + older_velocity_ * kDold;
 
@@ -184,6 +276,97 @@ Real3 CarTCell::CalculateDisplacement(const InteractionForce* force,
   return movement_at_next_step;
 }
 
+// Try to get attached to a tumor cell
+void CarTCell::TryToGetAttachedTo(TumorCell* victim, real_t squared_distance,
+                                  Random* rng) {
+  // If the tumor cell is not already attached to a CAR-T cell, is not dead and
+  // is not too far away.
+  if (!victim->IsAttachedToCart() && !victim->IsDead() &&
+      squared_distance < kSquaredMaxAdhesionDistanceCart) {
+    // factor of how high is the oncoprotein level of the cancer cell
+    real_t oncoprotein_scale_factor =
+        (victim->GetOncoproteinLevel() - kOncoproteinLimit) /
+        kOncoproteinDifference;
+    // Clamp scale_factor to be in [0,1]
+    if (oncoprotein_scale_factor > 1.0) {
+      oncoprotein_scale_factor = 1.0;
+    }
+    // If oncoprotein level is lower than the limit the cancer cell does not get
+    // detected
+    if (oncoprotein_scale_factor <= 0.0) {
+      // oncoprotein_scale_factor = 0.0; the probability is going to be 0 so
+      // return the function is the most efficient
+      return;
+    }
+
+    // factor of how far the cancer cell is
+    real_t distance_scale_factor =
+        (kMaxAdhesionDistanceCart - std::sqrt(squared_distance)) /
+        kDifferenceCartAdhesionDistances;
+    // Clamp scale_factor to be in [0,1]. We already checked that it is > 0
+    // because squared_distance < kSquaredMaxAdhesionDistanceCart
+    if (distance_scale_factor > 1.0) {
+      distance_scale_factor = 1.0;
+    }
+
+    // It tries to attach the CAR-T cell to the tumor cell with probability
+    // kAdhesionRateCart * oncoprotein_scale_factor * distance_scale_factor *
+    // kDtMechanics
+    if (rng->Uniform(0.0, 1.0) < kAdhesionRateCart * oncoprotein_scale_factor *
+                                     distance_scale_factor * kDtMechanics) {
+// avoid race condition. Only one cell can be attached to the tumor cell.
+#pragma omp critical
+      {
+        // We need to check again if the victim is not attached to a CAR-T cell
+        // yet. This could be made more efficiently with a semaphore for each
+        // cancer cell
+        if (!victim->IsAttachedToCart()) {
+          attached_to_tumor_cell_ = true;
+          attached_cell_ptr_ = victim->GetAgentPtr<TumorCell>();
+          victim->SetAttachedToCart(true);
+        }
+      }
+    }
+  }
+}
+
+// Try to induce apoptosis
+bool CarTCell::TryToInduceApoptosis(bdm::AgentPointer<TumorCell> attached_cell,
+                                    Random* rng) const {
+  // If there is no attached cell (this should not happen)
+
+  if (!attached_to_tumor_cell_) {
+    return false;
+  }
+
+  // factor of how high is the oncoprotein levelof the cancer cell
+  real_t scale_factor =
+      (attached_cell->GetOncoproteinLevel() - kOncoproteinLimit) /
+      kOncoproteinDifference;
+  // Clamp scale_factor to be in [0,1]
+  if (scale_factor > 1.0) {
+    scale_factor = 1.0;
+  }
+  // If oncoprotein level is lower than the limit the cancer cell does not
+  // become apoptotic
+  if (scale_factor < 0.0) {
+    // scale_factor = 0.0; the probability is going to be 0 so return the
+    // function is the most efficient
+    return false;
+  }
+  // CAR-T success of killing probability: aggressive cancer cells (high
+  // oncoprotein level) are more likely to be killed
+  const bool succeeded =
+      rng->Uniform(0.0, 1.0) < kKillRateCart * scale_factor * kDtMechanics;
+
+  // The CAR-T has succeeded to induce apoptosis on the Cancer Cell
+  if (succeeded) {
+    attached_cell->StartApoptosis();
+  }
+
+  return succeeded;
+}
+
 // Compute new oxygen or immunostimulatory factor concentration after
 // consumption/ secretion
 real_t CarTCell::ConsumeSecreteSubstance(int substance_id,
@@ -194,7 +377,7 @@ real_t CarTCell::ConsumeSecreteSubstance(int substance_id,
     res = (old_concentration + constant1_oxygen_) / constant2_oxygen_;
   } else if (substance_id ==
              immunostimulatory_factor_dgrid_->GetContinuumId()) {
-    // This point should never be reached
+    // CAR-T do not change immunostimulatory factor levels
     res = old_concentration;
   } else {
     throw std::invalid_argument("Unknown substance id: " +
@@ -230,7 +413,7 @@ void CarTCell::ComputeConstantsConsumptionSecretion() {
 
 /// Main behavior executed at each simulation step
 void StateControlCart::Run(Agent* agent) {
-  auto* sim = Simulation::GetActive();
+  Simulation* sim = Simulation::GetActive();
   // Run only every kDtCycle minutes, fmod does not work with the type
   // returned by GetSimulatedTime()
   if (sim->GetScheduler()->GetSimulatedSteps() % kStepsPerCycle != 0) {
@@ -241,7 +424,7 @@ void StateControlCart::Run(Agent* agent) {
     switch (cell->GetState()) {
       case CarTCellState::kAlive: {
         // the cell is growing to real_t its size before mitosis
-        // Probability of death= 1/CurrentLiveTime, division by 0
+        // Probability of death= 1/CurrentLiveTime,avoiding division by 0
         if (sim->GetRandom()->Uniform(1.0) <
             kDtCycle / std::max(cell->GetCurrentLiveTime(), kEpsilon)) {
           // the cell Dies
@@ -260,8 +443,8 @@ void StateControlCart::Run(Agent* agent) {
           cell->ComputeConstantsConsumptionSecretion();
           // Detach from tumor cell if it was attached
           if (cell->IsAttachedToTumorCell()) {
-            cell->GetAttachedCell()->SetAttachedToCart(false);
-            cell->SetAttachedCell(nullptr);
+            cell->GetAttachedCellPointer()->SetAttachedToCart(false);
+            cell->SetAttachedCellPointer(nullptr);
             cell->SetAttachedToTumorCell(false);
           }
         } else {
@@ -284,7 +467,7 @@ void StateControlCart::Run(Agent* agent) {
         // duration transition)
         if (kTimeApoptosis < cell->GetTimerState()) {
           // remove the cell from the simulation
-          auto* ctxt = sim->GetExecutionContext();
+          ExecutionContext* ctxt = sim->GetExecutionContext();
           ctxt->RemoveAgent(agent->GetUid());
         }
         break;