Multiple experiments and statistical analysis

Author: Lukas Breitwieser
In this tutorial we show how to collect and analyse data from multiple experiments.
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;

real_t gDivProb = 0.05;

We use the same simulation as in ST09-timeseries-plotting-basic. It is a simulation were agents divide with a specific division probability in each time step leading to exponential growth. We collect the number of agents in each time step. Have a look at ST09-timeseries-plotting-basic for more information.

We wrap the required simulation code in a function called Experiment which takes two parameters:

• the collected result data from a single invocation (output param)
• the division probability parameter

void Experiment(TimeSeries* result, real_t division_probability) {
    gDivProb = division_probability;
    
    auto set_param = [](Param* param) {
        param->simulation_time_step = 1.0;
    };
    Simulation simulation("MySimulation", set_param);
    
    StatelessBehavior rapid_division([](Agent* agent) {
      if (Simulation::GetActive()->GetRandom()->Uniform() < gDivProb) {
        bdm_static_cast<Cell*>(agent)->Divide();
      }
    });
    rapid_division.AlwaysCopyToNew();
    
    auto create_cell = [&](const Real3& position) {
      Cell* cell = new Cell(position);
      cell->SetDiameter(10);
      cell->AddBehavior(rapid_division.NewCopy());
      return cell;
    };
    
    simulation.GetResourceManager()->ClearAgents();
    ModelInitializer::CreateAgentsRandom(0, 200, 100, create_cell);
    simulation.GetScheduler()->FinalizeInitialization();
    
    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);
    
    simulation.GetScheduler()->Simulate(40);
    
    // move collected time series data from simulation to object result
    *result = std::move(*simulation.GetTimeSeries());
}

We want to run our experiment for 10 times with a different division probability parameter. We choose the division probability randomly between 0.04 and 0.06

std::vector<TimeSeries> individual_results(10);
Random rnd;
for(auto& ir : individual_results) {
    Experiment(&ir, rnd.Uniform(0.04, 0.06));
}

In the next step we want to combine the individual results. Therefore we calculate the mean, and min (error low), and max (error high) and store it in a merged TimeSeries object.

TimeSeries merged;
auto merger = [](const std::vector<real_t>& all_ys, 
                  real_t* y, real_t* eh, real_t* el) {
    *y = TMath::Mean(all_ys.begin(), all_ys.end());
    *el = *y - *TMath::LocMin(all_ys.begin(), all_ys.end());
    *eh = *TMath::LocMax(all_ys.begin(), all_ys.end()) - *y;
};
TimeSeries::Merge(&merged, individual_results, merger);

Now we can print the merged results and see how the simulations evolved over time, and how they differed from each other.

LineGraph g(&merged, "My result", "Time", "Number of agents", false, nullptr, 500, 300);
g.Add("num-agents", "Number of Agents", "LP", kBlue);
g.Draw();

Finally, let's fit an exponential function to the data.

g.GetTMultiGraph()->Fit("expo", "S");
g.Draw()

****************************************
Minimizer is Minuit2 / Migrad
Chi2                      =      0.85941
NDf                       =           37
Edm                       =  6.19237e-07
NCalls                    =           35
Constant                  =      4.67133   +/-   0.0158115   
Slope                     =    0.0503001   +/-   0.001853