Difference between revisions of "Team:Aix-Marseille/Collaborations"
(→Progamming code) |
(→Progamming code) |
||
| Line 48: | Line 48: | ||
|double V = 0.01; // Simulation volume with about 10^6 bacteria in ml | |double V = 0.01; // Simulation volume with about 10^6 bacteria in ml | ||
|double Vtot = 100.; // Volume of fermenter | |double Vtot = 100.; // Volume of fermenter | ||
| − | double D = 100.; // Dilution rate ml/hr | + | |double D = 100.; // Dilution rate ml/hr |
| − | double alpha = 3.4e-11;// growth yield g of substrate needed for 10^6 cells | + | |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 mumax = 3.0; // maximum growth rate on substrate doublings per hour |
| − | double KS = 0.1; // Monod constant for substrate (g/l) | + | |double KS = 0.1; // Monod constant for substrate (g/l) |
| − | + | | // | |
| − | + | | // Plasmid 1 parameters | |
| − | double KZ1 =100.0; // Growth inhibition constant | + | |double KZ1 =100.0; // Growth inhibition constant |
| − | int M1 = 1; // Hill constant for growth inhibition | + | |int M1 = 1; // Hill constant for growth inhibition |
| − | double k1 = 20.0; // Plasmid replication rate in hr^-1 | + | |double k1 = 20.0; // Plasmid replication rate in hr^-1 |
| − | double K1 = 0.1; // Plasmid replication inhibition constant | + | |double K1 = 0.1; // Plasmid replication inhibition constant |
| − | double Z1max = 10.0; // Maximum plsmid copy number | + | |double Z1max = 10.0; // Maximum plsmid copy number |
| − | + | | // | |
| − | + | | // Plasmid 2 parameters | |
| − | double KZ2 =100.0; // Growth inhibition constant | + | |double KZ2 =100.0; // Growth inhibition constant |
| − | int M2 = 1; // Hill constant for growth inhibition | + | |int M2 = 1; // Hill constant for growth inhibition |
| − | double k2 = 20.0; // Plasmid replication rate in hr^-1 | + | |double k2 = 20.0; // Plasmid replication rate in hr^-1 |
| − | double K2 = 0.1; // Plasmid replication inhibition constant | + | |double K2 = 0.1; // Plasmid replication inhibition constant |
| − | double Z2max = 10.0; // Maximum plsmid copy number | + | |double Z2max = 10.0; // Maximum plsmid copy number |
| − | + | | // | |
| − | double sigma = 5.0; // Division rate for large cells in hr^-1 | + | |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 ka1 = 2.0; // Toxicity parameter for toxin on plasmid 1 |
| − | double ka2 = 1.0; // Toxicity parameter for toxin on plasmid 2 | + | |double ka2 = 1.0; // Toxicity parameter for toxin on plasmid 2 |
| − | double kb = 2.0; // Ratio of anti-toxin to toxin production rates | + | |double kb = 2.0; // Ratio of anti-toxin to toxin production rates |
| − | + | | // | |
| − | + | | // Integrator parameters | |
| − | size_t total = 1e6; // Total number of bacteria at start | + | |size_t total = 1e6; // Total number of bacteria at start |
| − | double tmax = 100.0; // Number of hours to simulate | + | |double tmax = 100.0; // Number of hours to simulate |
| − | double dt = 0.05; // Timestep in hours | + | |double dt = 0.05; // Timestep in hours |
| − | double t = 0.0; // Current time | + | |double t = 0.0; // Current time |
| − | + | | | |
| − | double michaelis(double Vmax, double Km, double S) | + | |double michaelis(double Vmax, double Km, double S) |
| − | { | + | |{ |
| − | + | | return Vmax * S /(Km+S); | |
| − | } | + | |} |
| − | + | | | |
| − | gsl_rng *r; | + | |gsl_rng *r; |
| − | + | | | |
| − | void setup_seed() | + | |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; | + | |struct bstate { double Z[3]; } *population = NULL; |
| − | size_t pop_size = 0; | + | |size_t pop_size = 0; |
| − | + | | | |
| − | void pop_alloc(void) | + | |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) | + | |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) | + | |/* pop_delete(i) |
| − | + | | * overwrite population[i] with last element | |
| − | + | | * decrement pop_size | |
| − | + | | */ | |
| − | void pop_delete(size_t i) | + | |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[]) | + | |void printZ(double Z[]) |
| − | { | + | |{ |
| − | + | | for (int i = 0; i < 3; i++) | |
| − | + | | printf("%lf ", Z[i]); | |
| − | + | | printf("\n"); | |
| − | } | + | |} |
| − | + | | | |
| − | int main(void) | + | |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; | ||
}} | }} | ||
