Difference between revisions of "Team:Aix-Marseille/Collaborations"

(→‎Data recovery Bordeaux 2016)
(→‎Progamming code)
Line 21: Line 21:
  
 
===='''Progamming code''' ====
 
===='''Progamming code''' ====
 +
 +
{{hidden
 +
|Our computational model
 +
/*
 +
* Program to run a model of bacteria in a fermenter
 +
* with a 2 plasmid contention system.
 +
*/
 +
 +
#include <stdio.h>
 +
#include <stdlib.h>
 +
#include <string.h>
 +
#include <stdint.h>
 +
#include <sys/time.h>
 +
#include <gsl/gsl_rng.h>
 +
#include <gsl/gsl_randist.h>
 +
#include <gsl/gsl_sf.h>
 +
#include <math.h>
 +
 +
#define MAX_POPULATION 2E6
 +
 +
// Model parameters
 +
//
 +
// Fermentation variables
 +
double Sf    = 10.; // Substrate concentration (g/l) in feed solution
 +
double S      =  0.19; // Current substrate concentration in fementer (g/l)
 +
double V      =  0.01;  // Simulation volume with about 10^6 bacteria in ml
 +
double Vtot  = 100.; // Volume of fermenter
 +
double D      = 100.; // Dilution rate ml/hr
 +
double alpha  = 3.4e-11;// growth yield g of substrate needed for 10^6 cells
 +
double mumax  =  3.0; // maximum growth rate on substrate doublings per hour
 +
double KS    =  0.1; // Monod constant for substrate (g/l)
 +
    //
 +
    // Plasmid 1 parameters
 +
double KZ1    =100.0; // Growth inhibition constant
 +
int    M1    =    1; // Hill constant for growth inhibition
 +
double k1    = 20.0; // Plasmid replication rate in hr^-1
 +
double K1    =  0.1; // Plasmid replication inhibition constant
 +
double Z1max  = 10.0; // Maximum plsmid copy number
 +
    //
 +
    // Plasmid 2 parameters
 +
double KZ2    =100.0; // Growth inhibition constant
 +
int    M2    =    1; // Hill constant for growth inhibition
 +
double k2    = 20.0; // Plasmid replication rate in hr^-1
 +
double K2    =  0.1; // Plasmid replication inhibition constant
 +
double Z2max  = 10.0; // Maximum plsmid copy number
 +
    //
 +
double sigma  =  5.0; // Division rate for large cells in hr^-1
 +
    //
 +
    // Contention system parameters
 +
double ka1    =  2.0; // Toxicity parameter for toxin on plasmid 1
 +
double ka2    =  1.0; // Toxicity parameter for toxin on plasmid 2
 +
double kb    =  2.0; // Ratio of anti-toxin to toxin production rates
 +
    //
 +
    // Integrator parameters
 +
size_t total  = 1e6;  // Total number of bacteria at start
 +
double tmax  = 100.0; // Number of hours to simulate
 +
double dt    = 0.05; // Timestep in hours
 +
double t  = 0.0; // Current time
 +
 +
double michaelis(double Vmax, double Km, double S)
 +
{
 +
return Vmax * S /(Km+S);
 +
}
 +
 +
gsl_rng *r;
 +
 +
void setup_seed()
 +
{
 +
const gsl_rng_type * T;
 +
 +
gsl_rng_env_setup();
 +
T = gsl_rng_default;
 +
 +
struct timeval tv;
 +
gettimeofday(&tv,0);
 +
gsl_rng_default_seed = tv.tv_sec + tv.tv_usec;
 +
r = gsl_rng_alloc(T);
 +
}
 +
 +
struct bstate { double Z[3]; } *population = NULL;
 +
size_t pop_size = 0;
 +
 +
void pop_alloc(void)
 +
{
 +
population = calloc(MAX_POPULATION, sizeof(struct bstate));
 +
if (population == NULL) {
 +
fprintf(stderr, "%s: alloc error\n", __func__);
 +
exit(1);
 +
}
 +
}
 +
 +
void pop_append(struct bstate *p)
 +
{
 +
if (pop_size == MAX_POPULATION) {
 +
fprintf(stderr, "%s: max population reached\n", __func__);
 +
exit(1);
 +
}
 +
population[pop_size++] = *p;
 +
}
 +
 +
/* pop_delete(i)
 +
* overwrite population[i] with last element
 +
* decrement pop_size
 +
*/
 +
void pop_delete(size_t i)
 +
{
 +
if (i >= pop_size) {
 +
fprintf(stderr, "%s: out of range\n", __func__);
 +
exit(1);
 +
}
 +
 +
population[i] = population[--pop_size];
 +
}
 +
 +
void printZ(double Z[])
 +
{
 +
for (int i = 0; i < 3; i++)
 +
printf("%lf ", Z[i]);
 +
printf("\n");
 +
}
 +
 +
int main(void)
 +
{
 +
setup_seed();
 +
 +
/*
 +
* Create a random population of bacteria
 +
* in parameter space Z0(1-2),Z1(0-Z1max),Z2(0-Z2max)
 +
*/
 +
puts("Creating initial population");
 +
pop_alloc();
 +
 +
double AvZ[3] = {0};
 +
double growth = 0.0;
 +
uintmax_t free = 0;
 +
 +
for (size_t i = 0; i < total; i++) {
 +
struct bstate bacteria = {{
 +
1+gsl_rng_uniform(r), // Random size [1-2]
 +
Z1max, // Maximum number of plasmids
 +
Z2max // for both types
 +
}};
 +
pop_append(&bacteria);
 +
 +
AvZ[0] += bacteria.Z[0];
 +
AvZ[1] += bacteria.Z[1];
 +
AvZ[2] += bacteria.Z[2];
 +
}
 +
 +
puts("Starting integrator");
 +
 +
while (t < tmax) {
 +
/*
 +
* Display or output for visualization
 +
* population (density on Z1/Z2 (all Z0 or Z0>2),
 +
*/
 +
printf("%lf %lf %lf %zu %ju %lf %lf %lf\n",
 +
t, V, S, pop_size, free,
 +
AvZ[0]/pop_size, AvZ[1]/pop_size, AvZ[2]/pop_size);
 +
 +
memset(AvZ, 0, sizeof(double[3])); // average values for parameters
 +
growth = 0;
 +
free = 0; // number of bacteria with no plasmids
 +
 +
size_t i = pop_size;
 +
while (i-- != 0) {
 +
double* Z = population[i].Z;
 +
double cntrl  = gsl_rng_uniform(r) - (dt * D / Vtot);
 +
if (cntrl < 0) {
 +
pop_delete(i); // Bacterium washout
 +
continue;
 +
}
 +
// Bacteria grow
 +
AvZ[0] += Z[0];
 +
AvZ[1] += Z[1];
 +
AvZ[2] += Z[2];
 +
 +
if (Z[1] == 0.0 && Z[2] == 0.0)
 +
free++;
 +
 +
double dotZ0  = michaelis( mumax, KS, S ); // Growth on substrate
 +
dotZ0 *= michaelis( 1.0, gsl_sf_pow_int(Z[1], M1), KZ1 ); // Inhibition by plasmid 1
 +
dotZ0 *= michaelis( 1.0, gsl_sf_pow_int(Z[2], M2), KZ2 ); // Inhibition by plasmid 2
 +
double tox1 = gsl_sf_exp(-ka1*(Z[1]-kb*Z[2] ));
 +
double tox2 = gsl_sf_exp(-ka2*(Z[2]-kb*Z[1] ));
 +
dotZ0 *= fmin(1.0, tox1); // Inhibition by toxin 1
 +
dotZ0 *= fmin(1.0, tox2); // Inhibition by toxin 2
 +
 +
double dotZ1 = (Z[1] < 1.0)? 0.0 : michaelis( k1, K1, Z[0] ) * (Z1max - Z[1]);
 +
double dotZ2 = (Z[2] < 1.0)? 0.0 : michaelis( k2, K2, Z[0] ) * (Z2max - Z[2]);
 +
 +
Z[0] += dt * dotZ0; // Increment the internal state
 +
Z[1] += dt * dotZ1;
 +
Z[2] += dt * dotZ2;
 +
 +
growth += dt * dotZ0;
 +
 +
if (cntrl < (dt * sigma) && Z[0] > 2.0) {
 +
struct bstate new_bacterium = {{
 +
gsl_ran_gaussian_ziggurat(r, 0.05) + Z[0]/2.0,
 +
gsl_ran_binomial(r, 0.5, Z[1]),
 +
gsl_ran_binomial(r, 0.5, Z[2])
 +
}};
 +
for (int j = 0; j < 2; j++)
 +
Z[j] -= new_bacterium.Z[j];
 +
pop_append(&new_bacterium);
 +
}
 +
}
 +
 +
S += (D*dt*(Sf-S)/1000 - (growth*alpha*Vtot/V))/Vtot;
 +
 +
if (pop_size > 16e5) { // Too many bacteria
 +
// Throw out half of then and reduce the volume
 +
pop_size /= 2;
 +
V      /= 2.0;
 +
AvZ[0] /= 2.0;
 +
AvZ[1] /= 2.0;
 +
AvZ[2] /= 2.0;
 +
}
 +
if ((pop_size < 7e5) && (V < Vtot/2.0)) {
 +
// Too few bacteria increase the volume
 +
// and double the bacteria
 +
memcpy(population + pop_size, population, pop_size);
 +
pop_size *= 2;
 +
V      *= 2.0;
 +
AvZ[0] *= 2.0;
 +
AvZ[1] *= 2.0;
 +
AvZ[2] *= 2.0;
 +
}
 +
t += dt; // Increment time
 +
}
 +
// Output final population
 +
return 0;
 +
}
 +
 +
|style=text-align:center;
 +
}}
  
 
METTRE LE CODE DE FRANCOIS
 
METTRE LE CODE DE FRANCOIS

Revision as of 09:14, 17 October 2016