microbial-dynamics

# Microbial Dynamics: Population Dynamics & Modeling ## Overview Microbial Dynamics provides computational tools for modeling and analyzing microbial populations. This skill covers bacterial growth curve fitting using standard models (logistic, Gompertz, Baranyi), multi-species community dynamics via Lotka-Volterra equations, stochastic population simulation using the Gillespie algorithm, biofilm quantification from crystal violet assays, colony-forming unit enumeration with statistical analysis, bacterial genome annotation via Prokka, and simplified anaerobic digestion modeling. ## When to Use This Skill - Fitting bacterial growth curves from OD600 time-series data - Extracting growth parameters: lag phase duration, maximum growth rate (mu_max), carrying capacity (K) - Modeling multi-species microbial community interactions - Running stochastic simulations of gene expression or population dynamics - Processing crystal violet biofilm assay data - Calculating CFU/mL from serial dilution plating - Annotating bacterial genomes and extracting gene statistics - Modeling biogas production from anaerobic digestion **Related Skills:** For constraint-based metabolic modeling use `cobrapy`. For sequence manipulation and BLAST use `biopython`. For statistical analysis use `statistical-analysis`. ## Installation ```bash uv pip install scipy numpy pandas matplotlib ``` For genome annotation (optional): ```bash # conda install -c bioconda prokka ``` ## Quick Start ```python from scipy.optimize import curve_fit import numpy as np # Logistic growth model def logistic(t, y0, K, r, lag): return K / (1 + ((K - y0) / y0) * np.exp(-r * (t - lag))) # Example OD600 data time = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 16, 20, 24]) od600 = np.array([0.02, 0.02, 0.03, 0.06, 0.15, 0.38, 0.72, 1.05, 1.25, 1.42, 1.48, 1.50, 1.51, 1.51, 1.52]) popt, pcov = curve_fit(logistic, time, od600, p0=[0.02, 1.5, 0.5, 2.0], maxfev=10000) print(f"y0={popt[0]:.4f}, K={popt[1]:.3f}, r={popt[2]:.3f} h^-1, lag={popt[3]:.2f} h") ``` ## Core Capabilities ### 1. Growth Curve Modeling Fit OD600 data to standard microbial growth models. ```python import numpy as np from scipy.optimize import curve_fit def logistic(t, y0, K, r, lag): """Logistic growth model.""" return K / (1 + ((K - y0) / y0) * np.exp(-r * (t - lag))) def gompertz(t, y0, K, mu_max, lag): """Modified Gompertz growth model.""" return y0 + (K - y0) * np.exp(-np.exp((mu_max * np.e / (K - y0)) * (lag - t) + 1)) def baranyi(t, y0, K, mu_max, lag): """Baranyi growth model.""" A_t = t + (1 / mu_max) * np.log(np.exp(-mu_max * t) + np.exp(-mu_max * lag) - np.exp(-mu_max * (t + lag))) return K - np.log(1 + (np.exp(K) - np.exp(y0)) / np.exp(y0) * np.exp(-mu_max * A_t)) def fit_growth_curve(time, od600, model='logistic'): """Fit growth curve to OD data and return parameters.""" models = {'logistic': logistic, 'gompertz': gompertz, 'baranyi': baranyi} func = models[model] # Initial guesses y0_guess = od600[0] K_guess = od600.max() r_guess = 0.5 lag_guess = time[np.argmax(np.gradient(od600))] - 1 p0 = [y0_guess, K_guess, r_guess, max(lag_guess, 0)] bounds = ([0, 0, 0, 0], [np.inf, np.inf, 10, time.max()]) popt, pcov = curve_fit(func, time, od600, p0=p0, bounds=bounds, maxfev=10000) perr = np.sqrt(np.diag(pcov)) # R-squared residuals = od600 - func(time, *popt) ss_res = np.sum(residuals**2) ss_tot = np.sum((od600 - np.mean(od600))**2) r_squared = 1 - (ss_res / ss_tot) param_names = ['y0', 'K', 'mu_max', 'lag'] result = {name: {'value': val, 'std': err} for name, val, err in zip(param_names, popt, perr)} result['r_squared'] = r_squared result['model'] = model return result, popt # Compare models for model_name in ['logistic', 'gompertz', 'baranyi']: try: result, _ = fit_growth_curve(time, od600, model=model_name) print(f"{model_name}: R²={result['r_squared']:.4f}, " f"mu_max={result['mu_max']['value']:.3f} ± {result['mu_max']['std']:.3f} h⁻¹, " f"lag={result['lag']['value']:.2f} h") except RuntimeError: print(f"{model_name}: fitting failed") ``` ### 2. Lotka-Volterra Community Dynamics Simulate multi-species interactions. ```python import numpy as np from scipy.integrate import solve_ivp def lotka_volterra(t, N, r, K, alpha): """Generalized Lotka-Volterra for n species. Args: N: array of population sizes r: array of intrinsic growth rates K: array of carrying capacities alpha: interaction matrix (alpha[i,j] = effect of j on i) """ n = len(N) dNdt = np.zeros(n) for i in range(n): interaction = sum(alpha[i, j] * N[j] for j in range(n)) dNdt[i] = r[i] * N[i] * (1 - interaction / K[i]) return dNdt # 3-species community r = np.array([0.5, 0.4, 0.3]) # Growth rates K = np.array([1000, 800, 600]) # Carrying capacities alpha = np.array([ # Interaction matrix [1.0, 0.5, 0.1], # Species 1: self + competition from 2 and 3 [0.3, 1.0, 0.4], # Species 2 [0.2, 0.6, 1.0], # Species 3 ]) N0 = np.array([10, 10, 10]) # Initial populations sol = solve_ivp( lotka_volterra, [0, 100], N0, args=(r, K, alpha), t_eval=np.linspace(0, 100, 500), method='RK45' ) print("Final populations:") for i in range(3): print(f" Species {i+1}: {sol.y[i, -1]:.1f}") # Stability analysis: check if coexistence is stable A = np.diag(1/K) @ alpha try: eigenvalues = np.linalg.eigvals(A) stable = all(ev.real > 0 for ev in eigenvalues) print(f"Coexistence equilibrium is {'stable' if stable else 'unstable'}") except np.linalg.LinAlgError: print("Stability analysis failed") ``` ### 3. Stochastic Population Simulation Gillespie SSA for exact stochastic simulation. ```python import numpy as np def gillespie_ssa(propensity_func, stoich_matrix, x0, t_end, max_steps=100000): """Gillespie Stochastic Simulation Algorithm. Args: propensity_func: function(x) -> array of reaction propensities stoich_matrix: reactions x species stoichiometry matrix x0: initial state vector t_end: simulation end time """ t = 0 x = np.array(x0, dtype=float) times = [t] states = [x.copy()] for step in range(max_steps): props = propensity_func(x) total_prop = np.sum(props) if total_prop == 0 or t >= t_end: break # Time to next reaction dt = np.random.exponential(1 / total_prop) t += dt if t > t_end: break # Choose reaction reaction = np.searchsorted(np.cumsum(props), np.random.uniform(0, total_prop)) reaction = min(reaction, len(props) - 1) # Update state x += stoich_matrix[reaction] x = np.clip(x, 0, None) times.append(t) states.append(x.copy()) return np.array(times), np.array(states) # Birth-death-immigration model def propensities(x): N = x[0] birth_rate = 0.5 * N death_rate = 0.01 * N * (N - 1) # Density-dependent death immigration = 5.0 return np.array([birth_rate, death_rate, immigration]) stoich = np.array([[1], [-1], [1]]) # Birth: +1, Death: -1, Immigration: +1 # Run ensemble n_runs = 50 results = [] for i in range(n_runs): times, states = gillespie_ssa(propensities, stoich, [10], t_end=50) results.append((times, states)) # Summary statistics final_pops = [states[-1, 0] for _, states in results] print(f"Mean final population: {np.mean(final_pops):.1f} ± {np.std(final_pops):.1f}") ``` ### 4. Biofilm Quantification Process crystal violet biofilm assay data. ```python import numpy as np import pandas as pd def analyze_biofilm_cv(od_data, blank_od=0.05, conditions=None): """Analyze crystal violet biofilm assay. Args: od_data: dict of condition -> list of OD570 replicates blank_od: blank well OD for background subtraction """ results = [] for condition, replicates in od_data.items(): corrected = np.array(replicates) - blank_od corrected = np.clip(corrected, 0, None) results.append({ 'condition': condition, 'mean_od': np.mean(corrected), 'std_od': np.std(corrected, ddof=1), 'n': len(corrected), 'sem': np.std(corrected, ddof=1) / np.sqrt(len(corrected)) }) df = pd.DataFrame(results) # Normalize to control control_mean = df.iloc[0]['mean_od'] df['fold_change'] = df['mean_od'] / control_mean return df # Example: biofilm inhibition dose-response od_data = { 'Control': [0.85, 0.92, 0.88, 0.90], '1 uM': [0.80, 0.78, 0.82, 0.79], '10 uM': [0.55, 0.52, 0.58, 0.50], '100 uM': [0.20, 0.18, 0.22, 0.19], } df = analyze_biofilm_cv(od_data) print(df[['condition', 'mean_od', 'std_od', 'fold_change']]) ``` ### 5. CFU Enumeration Calculate colony-forming units from serial dilution plating. ```python import numpy as np from scipy import stats def calculate_cfu(counts, dilution_factor, volume_plated_ml=0.1): """Calculate CFU/mL from plate counts. Args: counts: list of colony counts per plate dilution_factor: dilution used (e.g., 1e-6 for 10^-6) volume_plated_ml: volume plated in mL """ counts = np.array(counts) cfu_per_ml = counts / (dilution_factor * volume_plated_ml) mean_cfu = np.mean(cfu_per_ml) std_cfu = np.std(cfu_per_ml, ddof=1) sem = std_cfu / np.sqrt(len(counts)) # 95% confidence interval ci = stats.t.interval(0.95, df=len(counts)-1, loc=mean_cfu, scale=sem) return { 'mean_cfu_per_ml': mean_cfu, 'std': std_cfu, 'sem': sem, 'ci_95': ci, 'n': len(counts), 'log10_cfu': np.log10(mean_cfu) } # Example result = calculate_cfu(counts=[42, 38, 45], dilution_factor=1e-6, volume_plated_ml=0.1) print(f"CFU/mL: {result['mean_cfu_per_ml']:.2e}") print(f"Log10 CFU/mL: {result['log10_cfu']:.2f}") print(f"95% CI: ({result['ci_95'][0]:.2e}, {result['ci_95'][1]:.2e})") ``` ### 6. Bacterial Genome Annotation Run Prokka and parse results. ```python import subprocess import pandas as pd def run_prokka(fasta_path, output_dir, prefix='genome', genus=None, species=None): """Annotate bacterial genome with Prokka.""" cmd = [ 'prokka', fasta_path,