Skip to content

Commit f3514d7

Browse files
vasvavvasvav
committed
Added 2 BioDynaMo-relevant examples in the "regulatory_networks" demo
1 parent ebb1910 commit f3514d7

3 files changed

Lines changed: 360 additions & 0 deletions

File tree

‎demo/regulatory_networks/src/bdm_ex1.h‎

Lines changed: 163 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,163 @@
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+
#ifndef BDM_EX1_H_
15+
#define BDM_EX1_H_
16+
17+
#include "biodynamo.h"
18+
19+
#include "core/behavior/regulatory_network.h"
20+
21+
namespace bdm {
22+
23+
enum Substances { kProtein };
24+
25+
#ifndef __ROOTCLING__
26+
struct ODE_system {
27+
const std::map<std::string, DiffusionGrid*>& mdg;
28+
const std::vector<real_t> param;
29+
30+
ODE_system(std::map<std::string, DiffusionGrid*>& m, const std::vector<real_t>& p)
31+
: mdg(m), param(p) {}
32+
33+
void operator()(const b_vector_t& x, b_vector_t& dxdt, real_t t,
34+
Agent* agent) const {
35+
auto& xyz = agent->GetPosition();
36+
auto dg = mdg.find("protein")->second;
37+
const real_t protein = dg->GetValue(xyz);
38+
39+
dxdt[0] = t * x[0] + param[0];
40+
dxdt[1] = t * x[1] * x[1] + param[1] * protein;
41+
dxdt[2] = t + x[2] + param[2];
42+
}
43+
};
44+
45+
struct ODE_jacobian {
46+
const std::map<std::string, DiffusionGrid*>& mdg;
47+
const std::vector<real_t> param;
48+
49+
ODE_jacobian(std::map<std::string, DiffusionGrid*>& m, const std::vector<real_t>& p)
50+
: mdg(m), param(p) {}
51+
52+
void operator()(const b_vector_t& x, b_matrix_t& jac, real_t t, b_vector_t& dfdt,
53+
Agent* agent) const {
54+
// auto& xyz = agent->GetPosition();
55+
// auto dg = mdg.find("protein")->second;
56+
// const real_t protein = dg->GetValue(xyz);
57+
58+
jac(0, 0) = t;
59+
jac(0, 1) = 0.0;
60+
jac(0, 2) = 0.0;
61+
jac(1, 0) = 0.0;
62+
jac(1, 1) = t * 2.0 * x[1];
63+
jac(1, 2) = 0.0;
64+
jac(2, 0) = 0.0;
65+
jac(2, 1) = 0.0;
66+
jac(2, 2) = 1.0;
67+
//
68+
dfdt[0] = dfdt[1] = dfdt[2] = 0.0;
69+
}
70+
};
71+
72+
struct ODE_output {
73+
void operator()(const b_vector_t& x, real_t t,
74+
const Agent* agent) {
75+
auto& xyz = agent->GetPosition();
76+
std::clog << agent->GetUid()
77+
<< ',' << xyz[0] << ',' << xyz[1] << ',' << xyz[2]
78+
<< ',' << t << ',' << x[0] << ',' << x[1] << ',' << x[2];
79+
std::clog << std::endl;
80+
}
81+
};
82+
#endif
83+
84+
namespace ex1 {
85+
86+
inline int Simulate(int argc, const char** argv) {
87+
std::ofstream fout("ex1.csv");
88+
// save the original buffer of std::clog
89+
std::streambuf* orig_clog_buff = std::clog.rdbuf();
90+
// redirect std::clog to point to the above file
91+
std::clog.rdbuf(fout.rdbuf());
92+
93+
// set-up the BioDynaMo simulation parameters
94+
auto set_parameters = [](Param* param) {
95+
param->output_dir = "ex1";
96+
param->use_progress_bar = false;
97+
param->bound_space = Param::BoundSpaceMode::kOpen;
98+
param->min_bound = 0.0;
99+
param->max_bound = +100.0;
100+
param->export_visualization = false;
101+
param->visualization_interval = 1;
102+
param->visualize_agents["Cell"] = { "diameter_", "volume_" };
103+
param->statistics = false;
104+
param->simulation_time_step = 1.0;
105+
param->visualize_diffusion = { Param::VisualizeDiffusion{"cytokine", true, true} };
106+
param->calculate_gradients = false;
107+
param->diffusion_method = "euler";
108+
};
109+
110+
Simulation sim(argc, argv, set_parameters);
111+
112+
auto* rm = sim.GetResourceManager();
113+
114+
// time-step of the BioDynaMo simulator
115+
const real_t dt_BDM = sim.GetParam()->simulation_time_step;
116+
// time-step of the regulatory network solver
117+
const real_t dt_RN = 0.01;
118+
// BioDynaMo's diffusion grid sample points in each dimension
119+
int n_DG = 51;
120+
121+
ModelInitializer::DefineSubstance(kProtein, "protein", 0.0/dt_BDM, 0.0/dt_BDM, n_DG);
122+
ModelInitializer::AddBoundaryConditions(
123+
kProtein, BoundaryConditionType::kNeumann,
124+
std::make_unique<ConstantBoundaryCondition>(0));
125+
126+
std::map<std::string, DiffusionGrid*> dg_map;
127+
dg_map.insert(std::make_pair("protein", rm->GetDiffusionGrid("protein")));
128+
129+
auto generate_cells = [&](const Real3& xyz) {
130+
Cell* c = new Cell();
131+
c->SetDiameter(1.0);
132+
c->SetAdherence(0.4);
133+
c->SetMass(1.0);
134+
c->SetPosition(xyz);
135+
#ifndef __ROOTCLING__
136+
c->AddBehavior(new RegulatoryNetwork(dt_RN, 1000, {1., 5., 7.},
137+
//ODE_solver::Euler,
138+
//ODE_solver::Rosenbrock,
139+
ODE_solver::RungeKutta,
140+
ODE_system(dg_map,{0.2,0.1,3.0}),
141+
ODE_jacobian(dg_map,{0.2,0.1,3.0}),
142+
ODE_output()));
143+
#endif
144+
return c;
145+
};
146+
147+
int n_cells = 1;
148+
ModelInitializer::CreateAgentsRandom(0.0, 0.0001, n_cells, generate_cells);
149+
150+
for (int s=0; s<10; s++)
151+
sim.GetScheduler()->Simulate(1);
152+
153+
// restore the original buffer of std::clog
154+
std::clog.rdbuf(orig_clog_buff);
155+
156+
return 0;
157+
}
158+
159+
} // namespace ex1
160+
161+
} // namespace bdm
162+
163+
#endif // BDM_EX1_H_

‎demo/regulatory_networks/src/bdm_ex2.h‎

Lines changed: 189 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,189 @@
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+
#ifndef BDM_EX2_H_
15+
#define BDM_EX2_H_
16+
17+
#include "biodynamo.h"
18+
19+
#include "core/behavior/regulatory_network.h"
20+
#include "core/agent/agent.h"
21+
22+
namespace bdm {
23+
24+
class MyCell : public Cell {
25+
BDM_AGENT_HEADER(MyCell, Cell, 1);
26+
27+
public:
28+
MyCell() {}
29+
explicit MyCell(const Real3& position) : Base(position), trail_(0.0) {}
30+
virtual ~MyCell() {}
31+
32+
void Initialize(const NewAgentEvent& event) override {
33+
Base::Initialize(event);
34+
35+
if (auto* mother = dynamic_cast<MyCell*>(event.existing_agent)) {
36+
if (event.GetUid() == CellDivisionEvent::kUid) {
37+
// copy properties from mother to daughter
38+
}
39+
}
40+
}
41+
42+
real_t GetTrail() const { return trail_; }
43+
void SetTrail(real_t t) { trail_ += t; }
44+
45+
private:
46+
/// keep track of the trail of the agent
47+
real_t trail_;
48+
};
49+
50+
#ifndef __ROOTCLING__
51+
struct Lorenz_rhs_ {
52+
void operator()(const b_vector_t& x, b_vector_t& dxdt, double t,
53+
Agent* agent) const {
54+
dxdt[0] = sigma * x[1] - sigma * x[0];
55+
dxdt[1] = rho * x[0] - x[1] - x[0] * x[2];
56+
dxdt[2] = -beta * x[2] + x[0] * x[1];
57+
}
58+
// https://en.wikipedia.org/wiki/Lorenz_system
59+
double sigma = 10.0;
60+
double rho = 28.0;
61+
double beta = 8.0 / 3.0;
62+
};
63+
64+
struct Lorenz_jac_ {
65+
void operator()(const b_vector_t& x, b_matrix_t& jac, double t, b_vector_t& dfdt,
66+
Agent* agent) const {
67+
jac(0, 0) = -sigma;
68+
jac(0, 1) = sigma;
69+
jac(0, 2) = 0.0;
70+
jac(1, 0) = rho - x[2];
71+
jac(1, 1) = -1.0;
72+
jac(1, 2) = -x[0];
73+
jac(2, 0) = x[1];
74+
jac(2, 1) = x[0];
75+
jac(2, 2) = -beta;
76+
//
77+
dfdt[0] = dfdt[1] = dfdt[2] = 0.0;
78+
}
79+
// https://en.wikipedia.org/wiki/Lorenz_system
80+
double sigma = 10.0;
81+
double rho = 28.0;
82+
double beta = 8.0 / 3.0;
83+
};
84+
85+
struct Lorenz_out_ {
86+
void operator()(const b_vector_t& x, real_t t, const
87+
Agent* agent) {
88+
auto& xyz = agent->GetPosition();
89+
std::clog << agent->GetUid()
90+
<< ',' << xyz[0] << ',' << xyz[1] << ',' << xyz[2]
91+
<< ',' << t << ',' << x[0] << ',' << x[1] << ',' << x[2];
92+
std::clog << std::endl;
93+
}
94+
};
95+
#endif
96+
97+
class Trajectory : public RegulatoryNetwork {
98+
BDM_BEHAVIOR_HEADER(Trajectory, RegulatoryNetwork, 1);
99+
100+
public:
101+
Trajectory() { AlwaysCopyToNew(); }
102+
#ifndef __ROOTCLING__
103+
Trajectory(real_t dt, int n_dt, const std::vector<real_t>& x)
104+
: RegulatoryNetwork(dt, n_dt, x, ODE_solver::Rosenbrock,
105+
Lorenz_rhs_(), Lorenz_jac_(), Lorenz_out_()) {}
106+
#endif
107+
virtual ~Trajectory() = default;
108+
109+
void Initialize(const NewAgentEvent& event) override {
110+
Base::Initialize(event);
111+
}
112+
113+
void Run(Agent* agent) override {
114+
#ifndef __ROOTCLING__
115+
Real3 xyz;
116+
for (int i=0; i<3; i++)
117+
xyz[i] = this->GetSpecie(i);
118+
119+
if (auto* cell = dynamic_cast<MyCell*>(agent)) {
120+
Real3 diff = xyz - cell->GetPosition();
121+
// calculate the cell trail
122+
cell->SetTrail(diff.Norm());
123+
// now set its new position
124+
cell->SetPosition(xyz);
125+
} else {
126+
Log::Fatal("Trajectory::Run", "agent is not of 'MyCell' type");
127+
}
128+
#endif
129+
}
130+
131+
};
132+
133+
namespace ex2 {
134+
135+
inline int Simulate(int argc, const char** argv) {
136+
std::ofstream fout("ex2.csv");
137+
// save the original buffer of std::clog
138+
std::streambuf* orig_clog_buff = std::clog.rdbuf();
139+
// redirect std::clog to point to the above file
140+
std::clog.rdbuf(fout.rdbuf());
141+
142+
// set-up the BioDynaMo simulation parameters
143+
auto set_parameters = [](Param* param) {
144+
param->output_dir = "ex2";
145+
param->use_progress_bar = false;
146+
param->bound_space = Param::BoundSpaceMode::kOpen;
147+
param->min_bound = 0.0;
148+
param->max_bound = +100.0;
149+
param->export_visualization = true;
150+
param->visualization_interval = 10;
151+
param->visualize_agents["MyCell"] = { "diameter_", "volume_", "trail_" };
152+
param->statistics = false;
153+
param->simulation_time_step = 1.0;
154+
};
155+
156+
Simulation sim(argc, argv, set_parameters);
157+
158+
// time-step of the regulatory network solver
159+
const real_t dt_RN = 0.01;
160+
161+
{
162+
Real3 xyz{1.0, 1.0, 1.0};
163+
164+
MyCell* c = new MyCell();
165+
c->SetDiameter(1.0);
166+
c->SetPosition(xyz);
167+
#ifndef __ROOTCLING__
168+
c->AddBehavior(new Trajectory(dt_RN, 10, {xyz[0], xyz[1], xyz[2]}));
169+
#else
170+
c->AddBehavior(new Trajectory());
171+
#endif
172+
173+
sim.GetExecutionContext()->AddAgent(c);
174+
}
175+
176+
for (int s=0; s<3000; s++)
177+
sim.GetScheduler()->Simulate(1);
178+
179+
// restore the original buffer of std::clog
180+
std::clog.rdbuf(orig_clog_buff);
181+
182+
return 0;
183+
}
184+
185+
} // namespace ex2
186+
187+
} // namespace bdm
188+
189+
#endif // BDM_EX2_H_

‎demo/regulatory_networks/src/regulatory_networks.h‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,9 @@
2424
#include "lorenz.h"
2525
#include "oscillator.h"
2626

27+
#include "bdm_ex1.h"
28+
#include "bdm_ex2.h"
29+
2730
namespace bdm {
2831

2932
inline int Simulate(int argc, const char** argv) {
@@ -35,6 +38,11 @@ inline int Simulate(int argc, const char** argv) {
3538
if (oscillator::Simulate(argc, argv))
3639
return 1;
3740

41+
if (bdm::ex1::Simulate(argc, argv))
42+
return 1;
43+
if (bdm::ex2::Simulate(argc, argv))
44+
return 1;
45+
3846
std::cout << "Simulation completed successfully!\n";
3947
return 0;
4048
}

0 commit comments

Comments
 (0)