The MATLAB code together with a representative output for the paper “Dispersal homogenizes communities via immigration even at low rates in a synthetic bacterial metacommunity” by Stilianos Fodelianakis, Alexander Lorz, Adriana ValenzuelaCuevas, Alan Barozzi, Jenny Booth and Daniele Daffonchio can be found here.
The folder consists of the following MATLAB functions:
 cells.m: represents the ordinary differential equation for the time dynamics of the strains as described in the Supplementary Material; it uses the growth rates computed in the function growth_r.m and the fraction of the vessel volume being transported per second, computed in the function frac_s.m.
 frac_s.m: computes the fraction of the vessel volume being transported per second at time t.
 growth_r.m: computes the growth rate of each of the three strains at time t.
 main.m: loops over all dispersal speeds and calculates for each dispersal speed the growthovermigration ratio averaged across strains, vessels and time.
 plottin.m: plots the results.
 sum_grow_over_immi.m: sums the growthoverimmigration ratio.
 temperatur.m: calculates the temperature in each of the three vessels as a function of time t.
In addition, there is the MATLAB data file KKX0.mat containing the length of the lag phase and the exponential growth rate, both fitted from the experimental data.
Running the MATLAB function main.m produces the following output:

 growth_over_immigration_speed_pump__3.25to3.85ul_per_s.csv: contains two rows; the first one with the dispersal speeds (in m^3/s) and the second one with the growthoverimmigration ratio.
 growth_over_immigration_speed_pump__3.25to3.85ul_per_s.jpg: shows the graph growthoverimmigration ratio as a function of dispersal speed.
 nnn_speed_pump__3.25to3.85ul_per_s.csv: contains 73 rows.Each column contains in this order:
 the dispersal speed,
 the population densities of strain B42 in vessel 25 at (30,60,90,120,150,180,210,240) minutes,
 the population densities of strain B42 in vessel 37 at (30,60,90,120,150,180,210,240) minutes,
 the population densities of strain B42 in vessel 42 at (30,60,90,120,150,180,210,240) minutes,
 the population densities of strain E310 in vessel 25 at (30,60,90,120,150,180,210,240) minutes,
 the population densities of strain E310 in vessel 37 at (30,60,90,120,150,180,210,240) minutes,
 the population densities of strain E310 in vessel 42 at (30,60,90,120,150,180,210,240) minutes,
 the population densities of strain E111 in vessel 25 at (30,60,90,120,150,180,210,240) minutes,
 the population densities of strain E111 in vessel 37 at (30,60,90,120,150,180,210,240) minutes,
 the population densities of strain E111 in vessel 42 at (30,60,90,120,150,180,210,240) minutes.
In each column, we have 1 dispersal speed, 3 strains, 3 vessels and 8 time points. This gives 1 + 3*3*8 = 73 rows.
The range of dispersal speeds in the first row of nnn_speed_pump__3.25to3.85ul_per_s.csv can be modified as follows.
In main.m, lines 71 – 84, we find the following code:
% large upper = speed_pump_arr(5); lower = speed_pump_arr(5)/25; stepsize=upper/50; % small upper = speed_pump_arr(5)/25; lower = 0; stepsize=upper/50; % medium upper = 3.85e9; lower = 3.25e9; stepsize=.3e9;
These are predefined ranges for the dispersal speeds consisting of upper bound, lower bound and step size. Commenting out the medium range selects the small range. Additionally, commenting out the small range selects the large range. Of course, a new range can also be defined, e.g. by overwriting the variables upper, lower, stepsize after the code defining the medium range in main.m.
Two additional scenarios can be tested by changing the code:
1) Scenario without heat coupling:
Change the file temperatur.m in the following way, while keeping everything else the same:
at the end of the script, replace:
dydt = f_tube.*[tem1y(1); tem2y(2); tem3y(3)]... %heat exchange +(t_delay<t)*(speed_pump*density/m)*([t_a;t_a;t_a]+... [ylag(3)t_a;ylag(1)t_a;ylag(2)t_a]*f_exp[y(1);y(2);y(3)]); %mixing
with
dydt = [0; 0; 0];
2) Scenario without immigration:
Change the file cells.m in the following way, while keeping everything else the same:
at the end of the script, replace
dndt = growth_rates'.*[n(1); n(2); n(3) ]... +(frac_speed)*[n(3)n(1);n(1)n(2);n(2)n(3)];
with
dndt = growth_rates'.*[n(1); n(2); n(3) ]%... % +(frac_speed)*[n(3)n(1);n(1)n(2);n(2)n(3)];