Simulation time series plotting and analysis

Author: Lukas Breitwieser
In this tutorial we show how to collect data during the simulation, and plot and analyse it at the end.
To this extent, we create a simulation where cells divide rapidly leading to exponential growth.

Let's start by setting up BioDynaMo notebooks.

%jsroot on
gROOT->LoadMacro("${BDMSYS}/etc/rootlogon.C");

INFO: Created simulation object 'simulation' with UniqueName='simulation'.

using namespace bdm::experimental;

auto set_param = [](Param* param) {
    param->simulation_time_step = 1.0;
};
Simulation simulation("MySimulation", set_param);

Let's create a behavior which divides cells with $5\%$ probability in each time step.
New cells should also get this behavior.
Therefore, we have to call AlwaysCopyToNew().
Otherwise, we would only see linear growth.

StatelessBehavior rapid_division([](Agent* agent) {
  if (Simulation::GetActive()->GetRandom()->Uniform() < 0.05) {
    bdm_static_cast<Cell*>(agent)->Divide();
  }
});
rapid_division.AlwaysCopyToNew();

Let's create a function that creates a cell at a specific position, with diameter = 10, and the rapid_division behavior.

auto create_cell = [](const Real3& position) {
  Cell* cell = new Cell(position);
  cell->SetDiameter(10);
  cell->AddBehavior(rapid_division.NewCopy());
  return cell;
};

As starting condition we want to create 100 cells randomly distributed in a cube with $min = 0, max = 200$

simulation.GetResourceManager()->ClearAgents();
ModelInitializer::CreateAgentsRandom(0, 200, 100, create_cell);
simulation.GetScheduler()->FinalizeInitialization();
VisualizeInNotebook();

Before we start the simulation, we have to tell BioDynaMo which data to collect.
We can do this with the TimeSeries::AddCollector function.
In this example we are interested in the number of agents ...

auto* ts = simulation.GetTimeSeries();
auto get_num_agents = [](Simulation* sim) {
  return static_cast<real_t>(sim->GetResourceManager()->GetNumAgents());
};
ts->AddCollector("num-agents", get_num_agents);

... and the number agents with $diameter < 5$.
We create a condition cond and pass it to an instance of Counter which calculates the number of agents for which cond(agent) evaluates to true.

auto cond = [](Agent* a) { return a->GetDiameter() < 5; };
ts->AddCollector("agents_lt_5", new bdm::experimental::Counter<real_t>(cond));

Now let's simulate 40 iterations

simulation.GetScheduler()->Simulate(40);

Now we can plot how the number of agents (in this case cells) and the number of agents with $diameter < 5$ evolved over time.

LineGraph g(ts, "my result", "Time", "Number of agents", 
                               true, nullptr, 500, 300);
g.Add("num-agents", "Number of Agents", "L", kBlue);
g.Add("agents_lt_5", "Number of Agents diam < 5", "L", kGreen);
g.Draw();

Let's try to fit an exponential function to verify our assumption that the cells grew exponentially.
Please visit the ROOT user-guide for more information regarding fitting

auto fitresult = g.GetTGraphs("num-agents")[0]->Fit("expo", "S");
g.Draw();

****************************************
Minimizer is Minuit2 / Migrad
Chi2                      =      377.117
NDf                       =           38
Edm                       =  2.80271e-11
NCalls                    =           49
Constant                  =      4.64791   +/-   0.00505256  
Slope                     =    0.0497729   +/-   0.000161219 

Indeed, the number of agents follow an exponential function $$y = \exp(slope * x + constant)$$ with constant = 4.6 and slope = 0.049 This corresponds to the division probability of $0.05$

This is how to change the color after the creation of g.
Also the position of the legend can be optimized.

g.GetTGraphs("num-agents")[0]->SetLineColor(kBlack);
g.SetLegendPos(1, 500, 20, 700);
g.Draw();

Let's save these results in multiple formats

g.SaveAs(Concat(simulation.GetOutputDir(), "/line-graph"), 
         {".root", ".svg", ".png", ".C"});