Difference between revisions of "Team:Aix-Marseille/Collaborations"
(→Progamming code) |
(→Progamming code) |
||
| Line 25: | Line 25: | ||
|Our computational model | |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; | |style=text-align:center; | ||
}} | }} | ||
