|
|
| Line 21: |
Line 21: |
| | | | |
| | ===='''Progamming code'''==== | | ===='''Progamming code'''==== |
| − | [[link=#code_collab|Show code]]
| |
| | | | |
| − | <html>
| + | Code de françois à ajouter |
| − | <style>
| + | |
| − | | + | |
| − | .code_collab{
| + | |
| − | display:none;
| + | |
| − | }
| + | |
| − | | + | |
| − | .code_collab:target{
| + | |
| − | display:block;
| + | |
| − | }
| + | |
| − | | + | |
| − | </style>
| + | |
| − | | + | |
| − | <div class="code_collab" id="code_collab">
| + | |
| − | <p>/*
| + | |
| − | * 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;
| + | |
| − | }</p></div>
| + | |
| − | </html>
| + | |
| | | | |
| | ===='''Results'''==== | | ===='''Results'''==== |
