gene_regulation.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// 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 CORE_BEHAVIOR_GENE_REGULATION_H_
#define CORE_BEHAVIOR_GENE_REGULATION_H_
 
#include <functional>
#include <vector>
 
#include "core/behavior/behavior.h"
#include "core/param/param.h"
#include "core/scheduler.h"
#include "core/simulation.h"
#include "core/util/root.h"
 
namespace bdm {
 
class GeneRegulation : public Behavior {
  BDM_BEHAVIOR_HEADER(GeneRegulation, Behavior, 1);
 
 public:
  GeneRegulation() { AlwaysCopyToNew(); }
 
  virtual ~GeneRegulation() = default;
 
  void Initialize(const NewAgentEvent& event) override {
    Base::Initialize(event);
 
    auto* other = event.existing_behavior;
    if (auto* gr = dynamic_cast<GeneRegulation*>(other)) {
      concentrations_ = gr->concentrations_;
      first_derivatives_ = gr->first_derivatives_;
    } else {
      Log::Fatal("GeneRegulation::EventConstructor",
                 "other was not of type GeneRegulation");
    }
  }
 
  void AddGene(const std::function<real_t(real_t, real_t)>& first_derivative,
               real_t initial_concentration) {
    first_derivatives_.push_back(first_derivative);
    concentrations_.push_back(initial_concentration);
  }
 
  const std::vector<real_t>& GetValues() const { return concentrations_; }
 
  void Run(Agent* agent) override {
    auto* sim = Simulation::GetActive();
    auto* param = sim->GetParam();
    auto* scheduler = sim->GetScheduler();
 
    const auto& timestep = param->simulation_time_step;
    uint64_t simulated_steps = scheduler->GetSimulatedSteps();
    const auto absolute_time = simulated_steps * timestep;
 
    if (param->numerical_ode_solver == Param::NumericalODESolver::kEuler) {
      // Euler
      for (uint64_t i = 0; i < first_derivatives_.size(); i++) {
        real_t slope = first_derivatives_[i](absolute_time, concentrations_[i]);
        concentrations_[i] += slope * timestep;
      }
    } else if (param->numerical_ode_solver == Param::NumericalODESolver::kRK4) {
      // Runge-Kutta 4
      for (uint64_t i = 0; i < first_derivatives_.size(); i++) {
        real_t interval_midpoint = absolute_time + timestep / 2.0;
        real_t interval_endpoint = absolute_time + timestep;
 
        real_t k1 = first_derivatives_[i](absolute_time, concentrations_[i]);
        real_t k2 = first_derivatives_[i](
            interval_midpoint, concentrations_[i] + timestep * k1 / 2.0);
        real_t k3 = first_derivatives_[i](
            interval_midpoint, concentrations_[i] + timestep * k2 / 2.0);
        real_t k4 = first_derivatives_[i](interval_endpoint,
                                          concentrations_[i] + timestep * k3);
 
        concentrations_[i] += timestep / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);
      }
    }
  }
 
 private:
  std::vector<real_t> concentrations_ = {};
 
  std::vector<std::function<real_t(real_t, real_t)>> first_derivatives_ = {};
};
 
}  // namespace bdm
 
#endif  // CORE_BEHAVIOR_GENE_REGULATION_H_