diff --git a/src/spent_media_model_interactions.py b/src/spent_media_model_interactions.py index b5a4ca1..4ea5d33 100644 --- a/src/spent_media_model_interactions.py +++ b/src/spent_media_model_interactions.py @@ -1,828 +1,824 @@ #!/usr/bin/env python # coding: utf-8 -# In[1]: +# In[2]: import sys import scipy import scipy.integrate as integr import numpy as np import pandas as pd import matplotlib as mpl import scipy.optimize as opt import matplotlib.pyplot as plt +from matplotlib import cm +from matplotlib.colors import ListedColormap, LinearSegmentedColormap -# In[2]: +# In[3]: ## Making sure the lag factor is working as expected t=list(range(1, 200)) lag_pow=2 comp=1 lag_factor0 = np.divide(np.power(t,lag_pow),np.add(np.power(t,lag_pow),np.power(0*comp,lag_pow))) lag_factor25 = np.divide(np.power(t,lag_pow),np.add(np.power(t,lag_pow),np.power(25*comp,lag_pow))) lag_factor50 = np.divide(np.power(t,lag_pow),np.add(np.power(t,lag_pow),np.power(50*comp,lag_pow))) lag_factor100 = np.divide(np.power(t,lag_pow),np.add(np.power(t,lag_pow),np.power(100*comp,lag_pow))) fig = plt.figure(figsize=[5,4], dpi=300) plt.plot(t,lag_factor0,lw=2,color='r') plt.plot(t,lag_factor25,lw=2,color='b') plt.plot(t,lag_factor50,lw=2,color='g') plt.plot(t,lag_factor100,lw=2,color='k') plt.xlabel("Time [a.u.]") plt.ylabel("Lag factor") plt.legend(('0','25','50','100'),loc='lower right') plt.show() outfile = '../figs/lag.pdf' fig.savefig(outfile,format='pdf',frameon='false',dpi=300,bbox_inches='tight') -# In[3]: +# In[4]: #### Based on scripts for the paper Piccardi et al. (2019) PNAS #### Available here https://gitlab.com/eccemic/facilitation2019 def f_monod(s,K,r): ###### Monod function f = rs/(s+K) # INPUT # s substrate concentration # K half-saturation of growth effect # r max growth rate f_num = np.multiply(s,r) f_den = np.add(s,K) ftemp = np.divide(f_num,f_den) # i.e. element-wise, f = rs/(s+K) return ftemp def f_Hill(s,K,r,h): ###### Hill function f = rs^h/(s^h+K^h) # INPUT # s substrate concentration # K half-saturation of growth effect # r max growth rate # h Hill coefficient spowh = np.power(s,h) Kpowh = np.power(K,h) f_num = np.multiply(spowh,r) # Numerator: r*s^h f_den = np.add(spowh,Kpowh) # Denominator: s^h+K^h ftemp = np.divide(f_num,f_den) return ftemp def f_mech(t,x,args): # Functional response, applies only to abundances in x[:-2] #### x is (n+4) length vector of abundances of n species ## 4 compound concentrations in last 4 positions r1max = args["r1"] # Max growth effect of compound 1 r2max = args["r2"] # Max growth effect of compound 2 r4max = args["r4"] # Max growth effect of compound 4 Kc1 = args["Kc1"] # Half-max effect of compound 1 Kc2 = args["Kc2"] # Half-max effect of compound 2 Kc4 = args["Kc4"] # Half-max effect of compound 4 Yinv1 = args["Y1"] # Compound 1 Yield [say gram of biomass per mol of nutrient] Yinv2 = args["Y2"] # Compound 2 Yield [say gram of biomass per mol of nutrient] Yinv4 = args["Y4"] # Compound 4 Yield [say gram of biomass per mol of nutrient] kprod3 = args["kprod3"] # Production rate of inhibitory compound 3 kprod4 = args["kprod4"] # Production rate of feeding compound 4 kupt = args["kup"] # Passive uptake rate of inhibitory compound 3 lagc1 = args["lc1"] # Length of lag phase when growing on compound 2 due to compound 1 lagc3 = args["lc3"] # Length of lag phase when growing on compound 2 due to compound 3 lagp = args["lp"] # Power factor for lag phase x_spec = x[:-4] c_comp1 = x[-4] c_comp2 = x[-3] c_comp3 = x[-2] c_comp4 = x[-1] t_eps=t+0.00001 # Prelims f_growth1 = f_monod(c_comp1,Kc1,r1max) f_growth2 = f_monod(c_comp2,Kc2,r2max) f_growth4 = f_monod(c_comp4,Kc4,r4max) f_lag_c1 = np.divide(np.power(t,lagp), np.add(np.power(t_eps,lagp),np.power(np.multiply(lagc1,c_comp1),lagp))) f_lag_c3 = np.divide(np.power(t,lagp), np.add(np.power(t_eps,lagp),np.power(np.multiply(lagc3,c_comp3),lagp))) Ns = len(x_spec) # Change in x due to growth through 3 nutrient compounds and one inhibitory compound # which delays growth dS_grow1 = np.multiply(f_growth1,x_spec) # Growth rc/(c+Kc)*x on comp 1 dS_grow2 = np.multiply(f_growth2,x_spec) # Growth rc/(c+Kc)*x on comp 2 dS_grow4 = np.multiply(f_growth4,x_spec) # Growth rc/(c+Kc)*x on comp 4 dS_lagc1 = f_lag_c1 dS_lagc3 = f_lag_c3 dC_comp1 = -1.0*np.dot(Yinv1, dS_grow1) dC_comp2 = -1.0*np.multiply(np.multiply(np.dot(Yinv2, dS_grow2), dS_lagc1), dS_lagc3) # lag only affects compound 2 uptake when # either compound 1 or compound 3 are present dC_comp3 = -1.0*np.dot(kupt,np.multiply(c_comp3,x_spec)) + np.dot(kprod3,x_spec) # passively taken up and/or produced dC_comp4 = np.add(-1.0*np.dot(Yinv4, dS_grow4), np.dot(kprod4,x_spec)) # produced and consumed # Assemble state vector dS = np.add(np.add(dS_grow1, np.multiply(np.multiply(dS_grow2, dS_lagc1),dS_lagc3)), dS_grow4) dX = np.append(np.append(np.append(np.append(dS,dC_comp1),dC_comp2), dC_comp3), dC_comp4) return dX def jac_df(s,K,r): # df_Monod(s)/ds = rK/(s+K)^2 df_numer = np.multiply(r,K) df_denom = np.power(np.add(s,K),2.0) return np.divide(df_numer,df_denom) def jac_df_Hill(s,K,r,h): # df_Hill(s)/ds = rK^h*hs^(h-1)/(s^h+K^h)^2 spowh = np.power(s,h) Kpowh = np.power(K,h) df_numer = h*np.multiply(r,Kpowh)*np.power(s,h-1.0) df_denom = np.power(np.add(spowh,Kpowh),2.0) return np.divide(df_numer,df_denom) def jac_df_lag(t,l,c,p): # df_lag(c)/dc = -p (tl)^p c^(p-1)/((lc)^p+t^p)^2 tpowp = np.power(t,p) lpowp = np.power(l,p) df_numer = -1.0*np.multiply(np.multiply(np.multiply(p,tpowp),lpowp),np.power(c,p-1)) df_denom = np.power(np.add(np.multiply(lpowp,np.power(c,p)),tpowp),2.0) return np.divide(df_numer,df_denom) def jac_mech(t, x, args): r1max = args["r1"] # Max growth effect of compound 1 r2max = args["r2"] # Max mortality of compound 2 r4max = args["r4"] # Max growth effect of compound 1 Kc1 = args["Kc1"] # Half-max effect of compound 1 Kc2 = args["Kc2"] # Half-max effect of compound 2 Kc4 = args["Kc4"] # Half-max effect of compound 1 Yinv1 = args["Y1"] # Compound Yield [say gram of biomass per mol of nutrient] Yinv2 = args["Y2"] # Compound Yield [say gram of biomass per mol of nutrient] Yinv4 = args["Y4"] # Compound Yield [say gram of biomass per mol of nutrient] kprod3 = args["kprod3"] # Production rate of inhibitory compound 3 kprod4 = args["kprod4"] # Production rate of feeding compound 4 kupt = args["kup"] # Passive uptake rate of inhibitory compound 3 lagc1 = args["lc1"] # Length of lag phase lagc3 = args["lc3"] # Length of lag phase lagp = args["lp"] # Power factor for lag phase x_spec = x[:-4] c_comp1 = x[-4] c_comp2 = x[-3] c_comp3 = x[-2] c_comp4 = x[-1] t_eps=t+0.00001 #### Prelims f_growth1 = f_monod(c_comp1,Kc1,r1max) f_growth2 = f_monod(c_comp2,Kc2,r2max) f_growth4 = f_monod(c_comp4,Kc4,r4max) f_lag_c1 = np.divide(np.power(t,lagp), np.add(np.power(t_eps,lagp),np.power(np.multiply(lagc1,c_comp1),lagp))) f_lag_c3 = np.divide(np.power(t,lagp), np.add(np.power(t_eps,lagp),np.power(np.multiply(lagc3,c_comp3),lagp))) Ns = len(x_spec) #### Compute partial derivatives wrt species abundance x and compoundconcs n1 and n2 # Effect on species dynamics: per-capita functional response # dfdS = np.subtract(f_growth,f_monod(c_tox,Kt,mmax)) dSdS = np.add(np.add(f_growth1, np.multiply(np.multiply(f_growth2,f_lag_c1),f_lag_c3), f_growth4)) dSdC1 = np.multiply(np.add(jac_df(c_comp1,Kc1,r1max)), np.multiply(np.multiply(f_growth2,jac_df_lag(t,lagc1,c_comp1,lagp)),f_lag_c3), x_spec) ## Check np.add term! dSdC2 = np.multiply(np.multiply(np.multiply(jac_df(c_comp2,Kc2,r2max),x_spec),f_lag_c1),f_lag_c3) dSdC3 = np.multiply(np.multiply(np.multiply(f_growth2,jac_df_lag(t,lagc3,c_comp3,lagp)),f_lag_c1),x_spec) dSdC4 = np.multiply(jac_df(c_comp4,Kc4,r4max),x_spec) # Effect on compound1 dynamics dC1dS = -1.0*np.multiply(Yinv1,f_growth1) # Compound 1 conc wrt species growth dC1dC1 = -1.0*np.dot(Yinv1,jac_df(c_comp1,Kc1,r1max)) # Compound conc 1 wrt nutrients dC1dC2 = 0.0 dC1dC3 = 0.0 dC1dC4 = 0.0 # Effect on compound 2 dynamics dC2dS = -1.0*np.multiply(np.multiply(np.multiply(Yinv2,f_growth2),f_lag_c1),f_lag_c3) dC2dC1 = -1.0*np.multiply(np.multiply(np.multiply(np.multiply(Yinv2,f_growth2),f_lag_c3), jac_df_lag(t,lagc1,c_comp1,lagp)),x_spec) dC2dC2 = -1.0*np.dot(Yinv2,dSdC2) dC2dC3 = -1.0*np.multiply(np.multiply(np.multiply(np.multiply(Yinv2,f_growth2),f_lag_c1), jac_df_lag(t,lagc3,c_comp3,lagp)),x_spec) dC2dC4 = 0.0 # Compound 3 does not change in concentration dC3dS = -kupt*c_comp3 + kprod3 dC3dC1 = 0.0 dC3dC2 = 0.0 dC3dC3 = -kupt*x_spec dC3dC4 = 0.0 # Effect on useful compound 4 dynamics dC4dS = np.add(-1.0*np.multiply(Yinv4,f_growth4),kprod4) dC4dC1 = 0.0 dC4dC2 = 0.0 dC4dC3 = 0.0 dC4dC4 = -1.0*np.dot(Yinv4,dSdC4) # Assemble Jacobian from partial derivatives dS = np.transpose(np.vstack((np.diag(dSdS),dSdC1,dSdC2,dSdC3,dSdC4))) dC1 = np.append(np.append(np.append(np.append(dC1dS,dC1dC1),dC1dC2),dC1dC3),dC1dC4) dC2 = np.append(np.append(np.append(np.append(dC2dS,dC2dC1),dC2dC2),dC2dC3),dC2dC4) dC3 = np.append(np.append(np.append(np.append(dC3dS,dC3dC1),dC3dC2),dC3dC3),dC3dC4) dC4 = np.append(np.append(np.append(np.append(dC4dS,dC4dC1),dC4dC2),dC4dC3),dC4dC4) J_out = np.vstack((dS,dC1,dC2,dC3,dC4)) return J_out def computeODE(N, t0, T, y0, pars_dict): Nspec = len(y0) x = integr.ode(f_mech,jac_mech).set_integrator('dopri5') xout = np.zeros([N+1,Nspec]) xout[0,:] = y0 x.set_initial_value(y0,t0).set_f_params(pars_dict).set_jac_params(pars_dict) dt = (T-t0)/N t = np.linspace(t0,T,N+1) idx_temp = 0 while x.successful() and idx_temp3: axtmp.plot(tplot,x_conc[:,3],ls='dashdot',lw=2,color='green') if legendtup: axtmp.legend(legendtup,loc='upper right') axtmp.set(xlabel="Time [a.u.]",ylabel="Abundance/conc. [a.u.]", title=title_text) return axtmp -# In[5]: +# In[6]: def draw_subplot_nolab(axtmp,tplot,xplot,legendtup,title_text,sp_color,num_C): x_spec = xplot[:,:-4] x_conc = xplot[:,-4:] axtmp.plot(tplot,x_spec,lw=2,color=sp_color) axtmp.plot(tplot,x_conc[:,0],ls='dashed',lw=2,color='k') axtmp.plot(tplot,x_conc[:,1],ls='dotted',lw=2,color='k') axtmp.plot(tplot,x_conc[:,2],ls='dashdot',lw=2,color='gray') if num_C>3: axtmp.plot(tplot,x_conc[:,3],ls='dashdot',lw=2,color='green') if legendtup: axtmp.legend(legendtup,loc='upper right') axtmp.set(xlabel="",ylabel="", title=title_text) return axtmp -# In[6]: +# In[101]: -def simulateGrowth(pars_s1, pars_s2, nc0, x0, C, t0, t_end): +def simulateGrowth(pars_s1, pars_s2, c0, nc0, x0, C, t0, t_end): #### Plot quartet tplot, xplots1 = computeODE(C, t0, t_end, np.hstack((x0[0],x0[1:])), pars_s1) x_conc_s1 = xplots1[:,-4:] tplot, xplots2_nc = computeODE(C, t0, t_end, np.hstack((x0[0],nc0)), pars_s2) tplot, xplots2 = computeODE(C, t0, t_end, np.hstack((x0[0],x0[1:])), pars_s2) tplot, xplots2_sp1 = computeODE(C, t0, t_end, np.hstack((x0[0],0.5*x_conc_s1[-1,:]+nc0)), pars_s2) tplot, xplots2_sp2 = computeODE(C, t0, t_end, np.hstack((x0[0],x_conc_s1[-1,:])), pars_s2) tplot, xplots2_sp3 = computeODE(C, t0, t_end, np.hstack((x0[0],0.5*x_conc_s1[-1,:]+c0)), pars_s2) #tplot, xplots2_sp4 = computeODE(C, t0, t_end, np.hstack((x0[0],0.5*x_conc_s1[-1,:]+c0-nc0)), pars_s2) return [tplot, xplots1, xplots2_nc, xplots2, xplots2_sp1, xplots2_sp2, xplots2_sp3]#, xplots2_sp4] -# In[7]: +# In[102]: def plotCurves(curves,num_C,outfile): #### Plot quartet fig, ax = plt.subplots(nrows=1, ncols=6, sharex='col', sharey='row', dpi=300, figsize=(12,2)) legendtup_s1 = (r'$S_1$',r'$C_1$',r'$C_2$',r'$C_3$',r'$C_4$') legendtup_s2 = (r'$S_2$',r'$C_1$',r'$C_2$',r'$C_3$',r'$C_4$') no_legend = "" ax[0].set_ylim(-1, 18) title_text = r"$S_1$ in MM" draw_subplot(ax[0], curves[0], curves[1], legendtup_s1, title_text, 'b',num_C) - title_text=r"$S_2$ in NC (A)" + title_text=r"$S_2$ in NC" draw_subplot(ax[1], curves[0], curves[2], legendtup_s2, title_text, 'r',num_C) - title_text=r"$S_2$ in MM (B)" + title_text=r"$S_2$ in MM" draw_subplot(ax[2], curves[0], curves[3], no_legend, title_text, 'r',num_C) - title_text=r"$S_2$ in 0.5" "\n" r"spent of $S_1$" "\n" r"+ NC (C)" + title_text=r"$S_2$ in 0.5" "\n" r"spent of $S_1$" "\n" r"+ NC" draw_subplot(ax[3], curves[0], curves[4], no_legend, title_text, 'r',num_C) - title_text=r"$S_2$ in" "\n" r"spent of $S_1$ (D)" + title_text=r"$S_2$ in" "\n" r"spent of $S_1$" draw_subplot(ax[4], curves[0], curves[5], no_legend, title_text, 'r',num_C) - title_text=r"$S_2$ in 0.5" "\n" r"spent of $S_1$" "\n" r"+ MM (E)" + title_text=r"$S_2$ in 0.5" "\n" r"spent of $S_1$" "\n" r"+ MM" draw_subplot(ax[5], curves[0], curves[6], no_legend, title_text, 'r',num_C) #title_text=r"$S_2$ in 0.5" "\n" r"spent of $S_1$" "\n" r"+ carbon" "\n" r"sources" #draw_subplot(ax[6], curves[0], curves[7], no_legend, title_text, 'r',num_C) plt.show() fig.savefig(outfile,format='pdf',dpi=300,bbox_inches='tight') -# In[8]: +# In[237]: -def plotAllCurves(curvesEC, curvesEC_IC, curvesEC_CF, curvesEC_CD, - curvesNS, curvesNS_IC, curvesNS_CF, curvesNS_CD, outfile): +def plotAllCurves(allcurves, outfile, c3_0): + + if c3_0>0: + curvesEC = allcurves[0] + curvesEC_IC = allcurves[1] + curvesEC_CF = allcurves[2] + curvesEC_CD = allcurves[3] + curvesNS = allcurves[4] + curvesNS_IC = allcurves[5] + curvesNS_CF = allcurves[6] + curvesNS_CD = allcurves[7] + else: + curvesEC = allcurves[0] + curvesEC_IC = allcurves[1] + curvesEC_CF = allcurves[2] + curvesEC_CD = allcurves[2] + curvesNS = allcurves[3] + curvesNS_IC = allcurves[4] + curvesNS_CF = allcurves[5] + curvesNS_CD = allcurves[5] + + #### Plot quartet fig, ax = plt.subplots(nrows=8, ncols=6, sharex='col', sharey='row', dpi=300, figsize=(12,15)) ax[0,0].set_ylim(-1, 18) num_C=3 legendtup_s1 = (r'$S_1$',r'$C_1$',r'$C_2$',r'$C_3$') legendtup_s2 = (r'$S_2$',r'$C_1$',r'$C_2$',r'$C_3$') no_legend = "" title_text = r"$S_1$ in MM" draw_subplot_nolab(ax[0,0], curvesEC[0], curvesEC[1], legendtup_s1, title_text, 'b', num_C) - title_text=r"$S_2$ in NC (A)" + title_text=r"$S_2$ in NC" draw_subplot_nolab(ax[0,1], curvesEC[0], curvesEC[2], legendtup_s2, title_text, 'r', num_C) - title_text=r"$S_2$ in MM (B)" + title_text=r"$S_2$ in MM" draw_subplot_nolab(ax[0,2], curvesEC[0], curvesEC[3], no_legend, title_text, 'r', num_C) - title_text=r"$S_2$ in 0.5" "\n" r"spent of $S_1$" "\n" r"+ NC (C)" + title_text=r"$S_2$ in 0.5" "\n" r"spent of $S_1$" "\n" r"+ NC" draw_subplot_nolab(ax[0,3], curvesEC[0], curvesEC[4], no_legend, title_text, 'r', num_C) - title_text=r"$S_2$ in" "\n" r"spent of $S_1$ (D)" + title_text=r"$S_2$ in" "\n" r"spent of $S_1$" draw_subplot_nolab(ax[0,4], curvesEC[0], curvesEC[5], no_legend, title_text, 'r', num_C) - title_text=r"$S_2$ in 0.5" "\n" r"spent of $S_1$" "\n" r"+ MM (E)" + title_text=r"$S_2$ in 0.5" "\n" r"spent of $S_1$" "\n" r"+ MM" draw_subplot_nolab(ax[0,5], curvesEC[0], curvesEC[6], no_legend, title_text, 'r', num_C) #title_text=r"$S_2$ in 0.5" "\n" r"spent of $S_1$" "\n" r"+ carbon" "\n" r"sources" #draw_subplot_nolab(ax[0,6], curvesEC[0], curvesEC[7], no_legend, title_text, 'r', num_C) title_text="" ax[1,0].set_ylim(-1, 18) legendtup_s1 = (r'$S_1$',r'$C_1$',r'$C_2$',r'$C_3$',r'$C_4$') legendtup_s2 = (r'$S_2$',r'$C_1$',r'$C_2$',r'$C_3$',r'$C_4$') draw_subplot_nolab(ax[1,0], curvesEC_IC[0], curvesEC_IC[1], legendtup_s1, title_text, 'b', num_C) draw_subplot_nolab(ax[1,1], curvesEC_IC[0], curvesEC_IC[2], legendtup_s2, title_text, 'r', num_C) draw_subplot_nolab(ax[1,2], curvesEC_IC[0], curvesEC_IC[3], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[1,3], curvesEC_IC[0], curvesEC_IC[4], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[1,4], curvesEC_IC[0], curvesEC_IC[5], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[1,5], curvesEC_IC[0], curvesEC_IC[6], no_legend, title_text, 'r', num_C) #draw_subplot_nolab(ax[1,6], curvesEC_IC[0], curvesEC_IC[7], no_legend, title_text, 'r', num_C) num_C=4 ax[2,0].set_ylim(-1, 18) legendtup_s1 = (r'$S_1$',r'$C_1$',r'$C_2$',r'$C_3$',r'$C_4$') legendtup_s2 = (r'$S_2$',r'$C_1$',r'$C_2$',r'$C_3$',r'$C_4$') draw_subplot_nolab(ax[2,0], curvesEC_CF[0], curvesEC_CF[1], legendtup_s1, title_text, 'b', num_C) draw_subplot_nolab(ax[2,1], curvesEC_CF[0], curvesEC_CF[2], legendtup_s2, title_text, 'r', num_C) draw_subplot_nolab(ax[2,2], curvesEC_CF[0], curvesEC_CF[3], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[2,3], curvesEC_CF[0], curvesEC_CF[4], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[2,4], curvesEC_CF[0], curvesEC_CF[5], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[2,5], curvesEC_CF[0], curvesEC_CF[6], no_legend, title_text, 'r', num_C) #draw_subplot_nolab(ax[1,6], curvesEC_IC[0], curvesEC_IC[7], no_legend, title_text, 'r', num_C) num_C=3 ax[3,0].set_ylim(-1, 18) legendtup_s1 = (r'$S_1$',r'$C_1$',r'$C_2$',r'$C_3$',r'$C_4$') legendtup_s2 = (r'$S_2$',r'$C_1$',r'$C_2$',r'$C_3$',r'$C_4$') draw_subplot_nolab(ax[3,0], curvesEC_CD[0], curvesEC_CD[1], legendtup_s1, title_text, 'b', num_C) draw_subplot_nolab(ax[3,1], curvesEC_CD[0], curvesEC_CD[2], legendtup_s2, title_text, 'r', num_C) draw_subplot_nolab(ax[3,2], curvesEC_CD[0], curvesEC_CD[3], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[3,3], curvesEC_CD[0], curvesEC_CD[4], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[3,4], curvesEC_CD[0], curvesEC_CD[5], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[3,5], curvesEC_CD[0], curvesEC_CD[6], no_legend, title_text, 'r', num_C) #draw_subplot_nolab(ax[1,6], curvesEC_IC[0], curvesEC_IC[7], no_legend, title_text, 'r', num_C) ax[4,0].set_ylim(-1, 18) legendtup_s1 = (r'$S_1$',r'$C_1$',r'$C_2$',r'$C_3$') legendtup_s2 = (r'$S_2$',r'$C_1$',r'$C_2$',r'$C_3$') draw_subplot_nolab(ax[4,0], curvesNS[0], curvesNS[1], legendtup_s1, title_text, 'b', num_C) draw_subplot_nolab(ax[4,1], curvesNS[0], curvesNS[2], legendtup_s2, title_text, 'r', num_C) draw_subplot_nolab(ax[4,2], curvesNS[0], curvesNS[3], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[4,3], curvesNS[0], curvesNS[4], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[4,4], curvesNS[0], curvesNS[5], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[4,5], curvesNS[0], curvesNS[6], no_legend, title_text, 'r', num_C) #draw_subplot_nolab(ax[2,6], curvesNS[0], curvesNS[7], no_legend, title_text, 'r', num_C) ax[5,0].set_ylim(-1, 18) legendtup_s1 = (r'$S_1$',r'$C_1$',r'$C_2$',r'$C_3$',r'$C_4$') legendtup_s2 = (r'$S_2$',r'$C_1$',r'$C_2$',r'$C_3$',r'$C_4$') draw_subplot_nolab(ax[5,0], curvesNS_IC[0], curvesNS_IC[1], legendtup_s1, title_text, 'b', num_C) draw_subplot_nolab(ax[5,1], curvesNS_IC[0], curvesNS_IC[2], legendtup_s2, title_text, 'r', num_C) draw_subplot_nolab(ax[5,2], curvesNS_IC[0], curvesNS_IC[3], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[5,3], curvesNS_IC[0], curvesNS_IC[4], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[5,4], curvesNS_IC[0], curvesNS_IC[5], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[5,5], curvesNS_IC[0], curvesNS_IC[6], no_legend, title_text, 'r', num_C) #draw_subplot_nolab(ax[5,6], curvesNS_CF[0], curvesNS_CF[7], no_legend, title_text, 'r', num_C) num_C=4 ax[6,0].set_ylim(-1, 18) legendtup_s1 = (r'$S_1$',r'$C_1$',r'$C_2$',r'$C_3$',r'$C_4$') legendtup_s2 = (r'$S_2$',r'$C_1$',r'$C_2$',r'$C_3$',r'$C_4$') draw_subplot_nolab(ax[6,0], curvesNS_CF[0], curvesNS_CF[1], legendtup_s1, title_text, 'b', num_C) draw_subplot_nolab(ax[6,1], curvesNS_CF[0], curvesNS_CF[2], legendtup_s2, title_text, 'r', num_C) draw_subplot_nolab(ax[6,2], curvesNS_CF[0], curvesNS_CF[3], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[6,3], curvesNS_CF[0], curvesNS_CF[4], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[6,4], curvesNS_CF[0], curvesNS_CF[5], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[6,5], curvesNS_CF[0], curvesNS_CF[6], no_legend, title_text, 'r', num_C) #draw_subplot_nolab(ax[3,6], curvesNS_CF[0], curvesNS_CF[7], no_legend, title_text, 'r', num_C) num_C=3 ax[7,0].set_ylim(-1, 18) draw_subplot_nolab(ax[7,0], curvesNS_CD[0], curvesNS_CD[1], legendtup_s1, title_text, 'b', num_C) draw_subplot_nolab(ax[7,1], curvesNS_CD[0], curvesNS_CD[2], legendtup_s2, title_text, 'r', num_C) draw_subplot_nolab(ax[7,2], curvesNS_CD[0], curvesNS_CD[3], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[7,3], curvesNS_CD[0], curvesNS_CD[4], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[7,4], curvesNS_CD[0], curvesNS_CD[5], no_legend, title_text, 'r', num_C) draw_subplot_nolab(ax[7,5], curvesNS_CD[0], curvesNS_CD[6], no_legend, title_text, 'r', num_C) #draw_subplot_nolab(ax[4,6], curvesNS_CD[0], curvesNS_CD[7], no_legend, title_text, 'r', num_C) plt.show() fig.savefig(outfile,format='pdf',dpi=300,bbox_inches='tight') -# In[9]: +# In[238]: def getRelativeAUC(curves): A=sum(curves[2][:,0]) # S2 in NC B=sum(curves[3][:,0]) # S2 in MM C=sum(curves[4][:,0]) # S2 in 0.5sp + NC D=sum(curves[5][:,0]) # S2 in sp E=sum(curves[6][:,0]) # S2 in 0.5sp + MM #F=sum(curves[7][:,0]) # S2 in 0.5sp + MM - NC return [(C-A)/(B-A), (D-A)/(B-A), (E-A)/(B-A)] -# In[10]: - - -#### Global settings -N = 500 -t0 = 0.0 -t_end = 500.0 # Days - -c3_0=2.0 # If environment is benign set to 0, otherwise to 2 -Nspec = 1 -x0s = [1.0] # Initial biomass -c0 = [1.0,1.0,c3_0,0.0] # Conc of compound 1 and 2 (nutrients) and 3 (inhibitory, taken up by sp1) -nc0 = [0.0,0.0,c3_0,0.0] # Non-carbonic medium -x0 = np.hstack((x0s,c0)) - -# Global constant -lp = 2 - -# Species 1 -r11 = 0.1 # Max growth effect of compound 1 -r12 = 0.0 # Max growth effect of compound 2 -r14 = 0.0 # Max growth effect of compound 4 -Kc11 = 1.0 # Half-max effect of compound 1 -Kc12 = 1.0 # Half-max effect of compound 2 -Kc14 = 1.0 # Half-max effect of compound 4 -Y11 = 0.1 # Compound 1 Yield [say gram of biomass per mol of nutrient] -Y12 = 0.1 # Compound 2 Yield [say gram of biomass per mol of nutrient] -Y14 = 0.1 # Compound 4 Yield [say gram of biomass per mol of nutrient] -lc11 = 0 # How much the concentration of compound 1 affects lag phase -lc13 = 0 # How much the concentration of compound 3 affects lag phase - -# Species 2 -Kc21 = Kc11 -Kc22 = Kc12 -Kc24 = Kc14 -Y21 = Y11 -Y22 = Y12 -Y24 = Y14 -lc21 = 0 # How much the concentration of nutrient 1 affects lag phase -kpr23 = 0.0 # Passive production of inhibitory compound 3 -kpr24 = 0.0 # Passive production of feeding compound 4 -kup23 = 0.0 # Passive uptake of compound 3 - -# Whether compound 3 prolongs the lag phase -lc23 = 200 # How much the concentration of compound 2 affects lag phase - - -# In[11]: +# In[239]: #### Exploitative competition +def get_curve_AUC(pars_s1, pars_s2, c0, nc0, x0, N, t0, t_end, num_C, filename): + curves=simulateGrowth(pars_s1,pars_s2,c0,nc0,x0,N,t0,t_end) + #plotCurves(curves,num_C,'../figs/'+filename) + return [curves, getRelativeAUC(curves)] -# Parameters for species 1 (first number) related to compounds (second number) -r11 = 0.0 # Max growth effect of compound 1 -r12 = 0.1 # Max growth effect of compound 2 -kpr13 = 0 # Passive production of compound 3 -kpr14 = 0 # Passive production of compound 4 -kup13 = 0 # Passive uptake of compound 3 - -pars_s1 = {'r1':r11,'r2':r12,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, - 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':kpr13,'kprod4':kpr14,'kup':kup13, - 'lc1':lc11,'lc3':lc13,'lp':lp} - -# Parameters for species 2 (first number) related to compounds (second number) -r21 = 0.0 # Max growth effect of compound 1 on species 1 -r22 = 0.1 # Max growth effect of compound 2 on species 1 -r24 = 0.0 # Max growth effect of compound 4 - -pars_s2 = {'r1':r21,'r2':r22,'r4':r24,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, - 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, - 'lc1':lc21,'lc3':lc23,'lp':lp} - -curvesEC=simulateGrowth(pars_s1,pars_s2,nc0,x0,N,t0,t_end) -plotCurves(curvesEC,3,'../figs/EC.pdf') -relativeCDE_EC = getRelativeAUC(curvesEC) - - -# In[12]: - - -#### Exploitative competition + interference competition - -# Parameters for species 1 (first number) related to compounds (second number) -r11 = 0.0 # Max growth effect of compound 1 -r12 = 0.1 # Max growth effect of compound 2 -kpr13 = 0.0005 # Passive production of compound 3 -kpr14 = 0 # Passive production of compound 4 -kup13 = 0 # Passive uptake of compound 3 - -pars_s1 = {'r1':r11,'r2':r12,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, - 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':kpr13,'kprod4':kpr14,'kup':kup13, - 'lc1':lc11,'lc3':lc13,'lp':lp} - -# Parameters for species 2 (first number) related to compounds (second number) -r21 = 0.0 # Max growth effect of compound 1 on species 1 -r22 = 0.1 # Max growth effect of compound 2 on species 1 -r24 = 0.0 # Max growth effect of compound 4 - -pars_s2 = {'r1':r21,'r2':r22,'r4':r24,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, - 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, - 'lc1':lc21,'lc3':lc23,'lp':lp} - -curvesEC_IC=simulateGrowth(pars_s1,pars_s2,nc0,x0,N,t0,t_end) -plotCurves(curvesEC_IC,4,'../figs/EC+IC.pdf') -relativeCDE_EC_IC = getRelativeAUC(curvesEC_IC) - - -# In[13]: +# In[240]: -#### Exploitative competition + crossfeeding -# Parameters for species 1 (first number) related to compounds (second number) -r11 = 0.0 -r12 = 0.1 -kpr13 = 0 # Passive production of compound 3 -kpr14 = 0.00005 # Passive production of compound 4 -kup13 = 0 # Passive uptake of compound 3 - -pars_s1 = {'r1':r11,'r2':r12,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, - 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':kpr13,'kprod4':kpr14,'kup':kup13, - 'lc1':lc11,'lc3':lc13,'lp':lp} - -# Parameters for species 2 (first number) related to compounds (second number) -r21 = 0.0 # Max growth effect of compound 1 -r22 = 0.1 # Max growth effect of compound 2 -r24 = 0.1 # Max growth effect of compound 4 - -pars_s2 = {'r1':r21,'r2':r22,'r4':r24,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, - 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, - 'lc1':lc21,'lc3':lc23,'lp':lp} - -curvesEC_CF=simulateGrowth(pars_s1,pars_s2,nc0,x0,N,t0,t_end) -plotCurves(curvesEC_CF,4,'../figs/EC+CF.pdf') -relativeCDE_EC_CF = getRelativeAUC(curvesEC_CF) - - -# In[14]: - - -#### Exploitative competition + crossdetoxification - -# Parameters for species 1 (first number) related to compounds (second number) -r11 = 0.0 -r12 = 0.1 -kpr13 = 0 # Passive production of compound 3 -kpr14 = 0 # Passive production of compound 4 -kup13 = 0.0001 # Passive uptake of compound 3 - -pars_s1 = {'r1':r11,'r2':r12,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, - 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':kpr13,'kprod4':kpr14,'kup':kup13, - 'lc1':lc11,'lc3':lc13,'lp':lp} - -# Parameters for species 2 (first number) related to compounds (second number) -r21 = 0.0 # Max growth effect of compound 1 -r22 = 0.1 # Max growth effect of compound 2 -r24 = 0 # Max growth effect of compound 4 - -pars_s2 = {'r1':r21,'r2':r22,'r4':r24,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, - 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, - 'lc1':lc21,'lc3':lc23,'lp':lp} - -curvesEC_CD=simulateGrowth(pars_s1,pars_s2,nc0,x0,N,t0,t_end) -plotCurves(curvesEC_CD,4,'../figs/EC+CD.pdf') -relativeCDE_EC_CD = getRelativeAUC(curvesEC_CD) - - -# In[15]: - - -#### Niche separation - -# Parameters for species 1 (first number) related to compounds (second number) -r11 = 0.1 # Max growth effect of compound 1 -r12 = 0.0 # Max growth effect of compound 2 -kpr13 = 0 # Passive production of compound 3 -kpr14 = 0 # Passive production of compound 4 -kup13 = 0 # Passive uptake of compound 3 - -pars_s1 = {'r1':r11,'r2':r12,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, - 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':kpr13,'kprod4':kpr14,'kup':kup13, - 'lc1':lc11,'lc3':lc13,'lp':lp} - -# Parameters for species 2 (first number) related to compounds (second number) -r21 = 0.0 # Max growth effect of compound 1 on species 1 -r22 = 0.1 # Max growth effect of compound 2 on species 1 -r24 = 0.0 # Max growth effect of compound 4 - -pars_s2 = {'r1':r21,'r2':r22,'r4':r24,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, - 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, - 'lc1':lc21,'lc3':lc23,'lp':lp} +#### Global settings +def runSimulations(kpr13_ic, kup13_cd, kpr14_cr, c3_0): +# kpr13_ic: If there is interference competition, this is how much of compound 3 is produced +# kup13_cd: If there is cross-detoxification, this is how much of compound 3 is taken up +# kpr14_cr: If there is cross-feeding, this is how much of compound 4 is produced +# c3_0: If environment is benign set to 0, otherwise to 2 + + N = 600 + t0 = 0.0 + t_end = 600.0 # Days + + Nspec = 1 + x0s = [1.0] # Initial biomass + c0 = [1.0,1.0,c3_0,0.0] # Conc of compound 1 and 2 (nutrients) and 3 (inhibitory, taken up by sp1) + nc0 = [0.0,0.0,c3_0,0.0] # Non-carbonic medium + x0 = np.hstack((x0s,c0)) + + # Global constant + lp = 2 + + # Species 1 + r11 = 0.1 # Max growth effect of compound 1 + r12 = 0.0 # Max growth effect of compound 2 + r14 = 0.0 # Max growth effect of compound 4 + Kc11 = 1.0 # Half-max effect of compound 1 + Kc12 = 1.0 # Half-max effect of compound 2 + Kc14 = 1.0 # Half-max effect of compound 4 + Y11 = 0.1 # Compound 1 Yield [say gram of biomass per mol of nutrient] + Y12 = 0.1 # Compound 2 Yield [say gram of biomass per mol of nutrient] + Y14 = 0.1 # Compound 4 Yield [say gram of biomass per mol of nutrient] + lc11 = 0 # How much the concentration of compound 1 affects lag phase + lc13 = 0 # How much the concentration of compound 3 affects lag phase + + # Species 2 + Kc21 = Kc11 + Kc22 = Kc12 + Kc24 = Kc14 + Y21 = Y11 + Y22 = Y12 + Y24 = Y14 + lc21 = 0 # How much the concentration of nutrient 1 affects lag phase + kpr23 = 0.0 # Passive production of inhibitory compound 3 + kpr24 = 0.0 # Passive production of feeding compound 4 + kup23 = 0.0 # Passive uptake of compound 3 + + + # Whether compound 3 prolongs the lag phase + lc23 = 200 # How much the concentration of compound 2 affects lag phase + + r12=0.1 + # Parameters for species 1 (first number) related to compounds (second number) + pars_s1 = {'r1':0,'r2':r12,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, + 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':0,'kprod4':0,'kup':0, + 'lc1':lc11,'lc3':lc13,'lp':lp} + # Parameters for species 2 (first number) related to compounds (second number) + pars_s2 = {'r1':0,'r2':0.1,'r4':0,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, + 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, + 'lc1':lc21,'lc3':lc23,'lp':lp} + [curves_EC, relCDE_EC ] = get_curve_AUC(pars_s1, pars_s2, c0, nc0, x0, N, t0, t_end, 3, '') + #----------------------------- + pars_s1 = {'r1':0,'r2':r12,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, + 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':kpr13_ic,'kprod4':0,'kup':0, + 'lc1':lc11,'lc3':lc13,'lp':lp} + pars_s2 = {'r1':0,'r2':0.1,'r4':0,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, + 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, + 'lc1':lc21,'lc3':lc23,'lp':lp} + [curves_EC_IC, relCDE_EC_IC] = get_curve_AUC(pars_s1, pars_s2, c0, nc0, x0, N, t0, t_end, 4, '') + #----------------------------- + pars_s1 = {'r1':0,'r2':r12,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, + 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':0,'kprod4':kpr14_cr,'kup':0, + 'lc1':lc11,'lc3':lc13,'lp':lp} + pars_s2 = {'r1':0,'r2':0.1,'r4':0.1,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, + 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, + 'lc1':lc21,'lc3':lc23,'lp':lp} + [curves_EC_CF, relCDE_EC_CF] = get_curve_AUC(pars_s1, pars_s2, c0, nc0, x0, N, t0, t_end, 4, '') + #----------------------------- + pars_s1 = {'r1':0,'r2':r12,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, + 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':0,'kprod4':0,'kup':kup13_cd, + 'lc1':lc11,'lc3':lc13,'lp':lp} + pars_s2 = {'r1':0,'r2':0.1,'r4':0,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, + 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, + 'lc1':lc21,'lc3':lc23,'lp':lp} + [curves_EC_CD, relCDE_EC_CD] = get_curve_AUC(pars_s1, pars_s2, c0, nc0, x0, N, t0, t_end, 4, '') + #----------------------------- + pars_s1 = {'r1':0.1,'r2':0,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, + 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':0,'kprod4':0,'kup':0, + 'lc1':lc11,'lc3':lc13,'lp':lp} + pars_s2 = {'r1':0,'r2':0.1,'r4':0,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, + 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, + 'lc1':lc21,'lc3':lc23,'lp':lp} + [curves_NS, relCDE_NS ] = get_curve_AUC(pars_s1, pars_s2, c0, nc0, x0, N, t0, t_end, 3, '') + #----------------------------- + pars_s1 = {'r1':0.1,'r2':0,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, + 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':kpr13_ic,'kprod4':0,'kup':0, + 'lc1':lc11,'lc3':lc13,'lp':lp} + pars_s2 = {'r1':0,'r2':0.1,'r4':0,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, + 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, + 'lc1':lc21,'lc3':lc23,'lp':lp} + [curves_NS_IC, relCDE_NS_IC] = get_curve_AUC(pars_s1, pars_s2, c0, nc0, x0, N, t0, t_end, 3, '') + + pars_s1 = {'r1':0.1,'r2':0,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, + 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':0,'kprod4':kpr14_cr,'kup':0, + 'lc1':lc11,'lc3':lc13,'lp':lp} + pars_s2 = {'r1':0,'r2':0.1,'r4':0.1,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, + 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, + 'lc1':lc21,'lc3':lc23,'lp':lp} + [curves_NS_CF, relCDE_NS_CF] = get_curve_AUC(pars_s1, pars_s2, c0, nc0, x0, N, t0, t_end, 4, '') + #----------------------------- + pars_s1 = {'r1':0.1,'r2':0,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, + 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':0,'kprod4':0,'kup':kup13_cd, + 'lc1':lc11,'lc3':lc13,'lp':lp} + pars_s2 = {'r1':0,'r2':0.1,'r4':0,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, + 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, + 'lc1':lc21,'lc3':lc23,'lp':lp} + [curves_NS_CD, relCDE_NS_CD] = get_curve_AUC(pars_s1, pars_s2, c0, nc0, x0, N, t0, t_end, 4, '') + + if c3_0>0: + all_curves = [curves_EC, curves_EC_IC, curves_EC_CF, curves_EC_CD, + curves_NS, curves_NS_IC, curves_NS_CF, curves_NS_CD] + all_relCDE = [relCDE_EC, relCDE_EC_IC, relCDE_EC_CF, relCDE_EC_CD, + relCDE_NS, relCDE_NS_IC, relCDE_NS_CF, relCDE_NS_CD] + else: + all_curves = [curves_EC, curves_EC_IC, curves_EC_CF, + curves_NS, curves_NS_IC, curves_NS_CF] + all_relCDE = [relCDE_EC, relCDE_EC_IC, relCDE_EC_CF, + relCDE_NS, relCDE_NS_IC, relCDE_NS_CF] + return[all_curves, all_relCDE] + + +# In[241]: + + +kpr13_ic = 0.0005 +kup13_cd = 0.00005 +kpr14_cr = 0.00001 +c3_0 = 0 #3.5 +[all_curves, AUC] = runSimulations(kpr13_ic, kup13_cd, kpr14_cr, c3_0) +plotAllCurves(all_curves, "../figs/20220609All.pdf", c3_0) + + +# In[228]: + + +# Plot AUCs with colored squares +from matplotlib import colors + +def colorFader(c1,c2,mix=0): #fade (linear interpolate) from color c1 (at mix=0) to c2 (mix=1) + c1=np.array(mpl.colors.to_rgb(c1)) + c2=np.array(mpl.colors.to_rgb(c2)) + return mpl.colors.to_hex((1-mix)*c1 + mix*c2) + +viridis = cm.get_cmap('viridis', 256) +newcolors = viridis(np.linspace(0, 1, 256)) + +c1='#d95f02' #orange +c2='#e6ab02' #yellow +c3='#1b9e77' #green +c4='#7570b3' #blue +n=64 +for x in range(n): + newcolors[x,:]=colors.to_rgba(colorFader(c1,c2,(x+1)/n)) + +m=192 +for x in range(m): + newcolors[n+x,:]=colors.to_rgba(colorFader(c3,c4,(x+1)/m)) + +newcmp = ListedColormap(newcolors, name='two_gradients') -curvesNS=simulateGrowth(pars_s1,pars_s2,nc0,x0,N,t0,t_end) -plotCurves(curvesNS,3,'../figs/NS.pdf') -relativeCDE_NS = getRelativeAUC(curvesNS) +kpr13_ic = 0.0005 +kup13_cd = 0.00005 +kpr14_cr = 0.00001 +c3_0 = 3.5 +[all_curves, AUC] = runSimulations(kpr13_ic, kup13_cd, kpr14_cr, c3_0) +fig, axs = plt.subplots(1, 1, figsize=(1.7, 3/8*(6.5+c3_0*0.75)), constrained_layout=True, squeeze=False) -# In[16]: +axs[0,0].pcolormesh(AUC, cmap=newcmp, rasterized=True, vmin=0, vmax=4, edgecolors='black') +axs[0,0].invert_yaxis() +fig.colorbar(psm, ax=axs[0,0]) +fig.savefig('../figs/harsh_model.pdf',format='pdf',dpi=300,bbox_inches='tight') -#### Niche separation + interference competition +AUC -# Parameters for species 1 (first number) related to compounds (second number) -r11 = 0.1 # Max growth effect of compound 1 -r12 = 0.0 # Max growth effect of compound 2 -kpr13 = 0.0005 # Passive production of compound 3 -kpr14 = 0 # Passive production of compound 4 -kup13 = 0 # Passive uptake of compound 3 -pars_s1 = {'r1':r11,'r2':r12,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, - 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':kpr13,'kprod4':kpr14,'kup':kup13, - 'lc1':lc11,'lc3':lc13,'lp':lp} +# In[151]: -# Parameters for species 2 (first number) related to compounds (second number) -r21 = 0.0 # Max growth effect of compound 1 on species 1 -r22 = 0.1 # Max growth effect of compound 2 on species 1 -r24 = 0.0 # Max growth effect of compound 4 -pars_s2 = {'r1':r21,'r2':r22,'r4':r24,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, - 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, - 'lc1':lc21,'lc3':lc23,'lp':lp} +# Plot experimental data +df = pd.read_excel(r'../AUC_OD_NewNorm.xlsx', sheet_name='Summary') +df_At = df[df["GrowingSpecies"] == "At"][["SM/2+NC", "SM", "SM/2+MM"]] +df_Ct = df[df["GrowingSpecies"] == "Ct"][["SM/2+NC", "SM", "SM/2+MM"]] +df_Oa = df[df["GrowingSpecies"] == "Oa"][["SM/2+NC", "SM", "SM/2+MM"]] -curvesNS_IC=simulateGrowth(pars_s1,pars_s2,nc0,x0,N,t0,t_end) -plotCurves(curvesNS_IC,3,'../figs/NS+IC.pdf') -relativeCDE_NS_IC = getRelativeAUC(curvesNS_IC) +fig, axs = plt.subplots(3, 1, figsize=(1.4, 5), constrained_layout=True, squeeze=False) +axs[0,0].pcolormesh(df_At.to_numpy(), cmap=newcmp, rasterized=True, vmin=0, vmax=4, edgecolors='black') +axs[0,0].invert_yaxis() -# In[17]: +axs[1,0].pcolormesh(df_Ct.to_numpy(), cmap=newcmp, rasterized=True, vmin=0, vmax=4, edgecolors='black') +axs[1,0].invert_yaxis() +print(df_Ct) +axs[2,0].pcolormesh(df_Oa.to_numpy(), cmap=newcmp, rasterized=True, vmin=0, vmax=4, edgecolors='black') +axs[2,0].invert_yaxis() -#### Niche separation + crossfeeding +fig.savefig('../figs/At_Ct_Oa_AUCs.pdf',format='pdf',dpi=300,bbox_inches='tight') -# Parameters for species 1 (first number) related to compounds (second number) -kpr13 = 0 # Passive production of compound 3 -kpr14 = 0.00005 # Passive production of compound 4 -kup13 = 0 # Passive uptake of compound 3 -pars_s1 = {'r1':r11,'r2':r12,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, - 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':kpr13,'kprod4':kpr14,'kup':kup13, - 'lc1':lc11,'lc3':lc13,'lp':lp} +# In[224]: -# Parameters for species 2 (first number) related to compounds (second number) -r21 = 0.0 # Max growth effect of compound 1 -r22 = 0.1 # Max growth effect of compound 2 -r24 = 0.1 # Max growth effect of compound 4 -pars_s2 = {'r1':r21,'r2':r22,'r4':r24,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, - 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, - 'lc1':lc21,'lc3':lc23,'lp':lp} +# Parameter sweep for kpr14_cr (production of cross-fed compound) +kpr13_ic = 0.0005 +kup13_cd = 0.00005 +kpr14_cr = [0.000005, 0.00001, 0.00005, 0.0001, 0.0002] +c3_0 = 3.5 -curvesNS_CF=simulateGrowth(pars_s1,pars_s2,nc0,x0,N,t0,t_end) -plotCurves(curvesNS_CF,4,'../figs/NS+CF.pdf') -relativeCDE_NS_CF = getRelativeAUC(curvesNS_CF) +fig, axs = plt.subplots(1, len(kpr14_cr), figsize=(1.5 * len(kpr14_cr), 3), + constrained_layout=True, squeeze=False) +for i in range(0, len(kpr14_cr)): + [all_curves, AUC] = runSimulations(kpr13_ic, kup13_cd, kpr14_cr[i], c3_0) -# In[18]: + psm = axs[0,i].pcolormesh(AUC, cmap=newcmp, rasterized=True, vmin=0, vmax=4, edgecolors='black') + axs[0,i].invert_yaxis() + axs[0,i].title.set_text('p_1,4='+str(kpr14_cr[i])) +fig.colorbar(psm, ax=axs[0,len(kpr14_cr)-1]) +fig.savefig('../figs/sweep_kpr14/sweep_kpr14_comparison.pdf',format='pdf',dpi=300,bbox_inches='tight') -#### Niche separation + crossdetoxification -# Parameters for species 1 (first number) related to compounds (second number) -r11 = 0.1 # Max growth effect of compound 1 -r12 = 0.0 # Max growth effect of compound 2 -kpr13 = 0 # Passive production of compound 3 -kpr14 = 0 # Passive production of compound 4 -kup13 = 0.0001 # Passive uptake of compound 3 +# In[225]: -pars_s1 = {'r1':r11,'r2':r12,'r4':r14,'Kc1':Kc11,'Kc2':Kc12,'Kc4':Kc14, - 'Y1':Y11,'Y2':Y12,'Y4':Y14,'kprod3':kpr13,'kprod4':kpr14,'kup':kup13, - 'lc1':lc11,'lc3':lc13,'lp':lp} -# Parameters for species 2 (first number) related to compounds (second number) -r21 = 0.0 # Max growth effect of compound 1 -r22 = 0.1 # Max growth effect of compound 2 -r24 = 0 # Max growth effect of compound 4 +# Parameter sweep for kpr13_cr (production of inhibitory compound) +kpr13_ic = [0.000125, 0.00025, 0.0005, 0.001, 0.0025] +kup13_cd = 0.00005 +kpr14_cr = 0.00001 +c3_0 = 3.5 -pars_s2 = {'r1':r21,'r2':r22,'r4':r24,'Kc1':Kc21,'Kc2':Kc22,'Kc4':Kc24, - 'Y1':Y21,'Y2':Y22,'Y4':Y24,'kprod3':kpr23,'kprod4':kpr24,'kup':kup23, - 'lc1':lc21,'lc3':lc23,'lp':lp} +fig, axs = plt.subplots(1, len(kpr13_ic), figsize=(1.5 * len(kpr13_ic), 3), + constrained_layout=True, squeeze=False) -curvesNS_CD=simulateGrowth(pars_s1,pars_s2,nc0,x0,N,t0,t_end) -plotCurves(curvesNS_CD,4,'../figs/NS+CD.pdf') -relativeCDE_NS_CD = getRelativeAUC(curvesNS_CD) +for i in range(0, len(kpr13_ic)): + [all_curves, AUC] = runSimulations(kpr13_ic[i], kup13_cd, kpr14_cr, c3_0) + psm = axs[0,i].pcolormesh(AUC, cmap=newcmp, rasterized=True, vmin=0, vmax=4, edgecolors='black') + axs[0,i].invert_yaxis() + axs[0,i].title.set_text('p_1,3='+str(kpr13_ic[i])) -# In[19]: +fig.colorbar(psm, ax=axs[0,len(kpr13_ic)-1]) +fig.savefig('../figs/sweep_kpr13_comparison.pdf',format='pdf',dpi=300,bbox_inches='tight') -# AUC of all interactions +# In[229]: -if c3_0==0: - plotAllCurves(curvesEC, curvesEC_IC, curvesEC_CF, curvesEC_CD, - curvesNS, curvesNS_IC, curvesNS_CF, curvesNS_CD, '../figs/All.pdf') - AUC=[relativeCDE_EC, relativeCDE_EC_IC, relativeCDE_EC_CF, relativeCDE_EC_CD, - relativeCDE_NS, relativeCDE_NS_IC, relativeCDE_NS_CF, relativeCDE_NS_CD] - print(AUC) -else: - plotAllCurves(curvesEC, curvesEC_IC, curvesEC_CF, curvesEC_CD, - curvesNS, curvesNS_IC, curvesNS_CF, curvesNS_CD, '../figs/All_inhib.pdf') - AUC_inhib=[relativeCDE_EC, relativeCDE_EC_IC, relativeCDE_EC_CF, relativeCDE_EC_CD, - relativeCDE_NS, relativeCDE_NS_IC, relativeCDE_NS_CF, relativeCDE_NS_CD] - print(AUC_inhib) +# Parameter sweep for kup13_cr (uptake of inhibitory compound) +kpr13_ic = 0.0005 +kup13_cd = [0.000005, 0.00001, 0.00005, 0.0001, 0.0002] +kpr14_cr = 0.00001 +c3_0 = 3.5 -# In[20]: +fig, axs = plt.subplots(1, len(kup13_cd), figsize=(1.5 * len(kup13_cd), 3), + constrained_layout=True, squeeze=False) +for i in range(0, len(kup13_cd)): + [all_curves, AUC] = runSimulations(kpr13_ic, kup13_cd[i], kpr14_cr, c3_0) -## Non-inhibitory environment -# Interaction=[C, D, E] -# EC=[0.0, 0.0, 1.0] -# EC+IC=[0.0, 0.0, 0.6148869705931375] -# EC+CF=[0.11943890903962877, 0.245419510551685, 1.1388376113005] -# NS=[0.4924147048163176, 1.0, 1.5074901036674009] -# NS+IC=[0.29410504357843814, 0.3708855372591011, 0.9351830736554286] -# NS+CF=[0.6274933023154798, 1.2758809375070677, 1.6490335573759083] + psm = axs[0,i].pcolormesh(AUC, cmap=newcmp, rasterized=True, vmin=0, vmax=4, edgecolors='black') + axs[0,i].invert_yaxis() + axs[0,i].title.set_text('u_1,3='+str(kup13_cd[i])) -## Inhibitory environment -# EC=[0.0, 0.0, 0.6298844610312351] -# EC+IC=[0.0, 0.0, 0.2862847338927636] -# EC+CF=[0.25450859644725093, 0.5229566786824308, 1.1148227458307745] -# EC+CD=[0.0, 0.0, 0.7710933686130353] -# NS=[0.2821971649816627, 1.0, 0.9826016271674324] -# NS+IC=[0.13943237908543293, 0.24192509125585673, 0.42807890918642505] -# NS+CF=[0.6723586130218301, 1.8106816153523129, 1.5479069980498832] -# NS+CD=[0.34955458756258073, 1.3441574243064927, 1.1936229033557797] +fig.colorbar(psm, ax=axs[0,len(kup13_cd)-1]) +fig.savefig('../figs/sweep_kup13_comparison.pdf',format='pdf',dpi=300,bbox_inches='tight') -# In[21]: +# In[187]: def plotSpecies2(curves, legendtup, title_text, outfile): #### Plot single panel with all curves of species 2 fig = plt.figure(dpi=300, figsize=(3,3)) ax = fig.add_subplot(1, 1, 1) ax.set_ylim(-1, 18) ax.plot(curves[0],curves[2][:,:-4],lw=2,color='gray') ax.plot(curves[0],curves[3][:,:-4],lw=2,color='k') ax.plot(curves[0],curves[4][:,:-4],ls='dotted',lw=2,color='#a6761d') ax.plot(curves[0],curves[5][:,:-4],lw=2,color='#a6761d') ax.plot(curves[0],curves[6][:,:-4],ls='dashed',lw=2,color='#a6761d') ax.legend(legendtup,loc='upper left') ax.set(xlabel="Time [a.u.]",ylabel="Abundance/conc. [a.u.]", title=title_text) plt.show() fig.savefig(outfile,format='pdf',dpi=300,bbox_inches='tight') -# In[22]: +# In[223]: + +# Generate figure to compare curves under different scenarios +curvesNS_CF = all_curves[6] +curvesNS_CD = all_curves[7] -legendtup = (r'A',r'B',r'C',r'D',r'E') +legendtup = (r'NC',r'MM',r'SM/2+NC',r'SM',r'SM/2+MM') plotSpecies2(curvesNS_CF, legendtup, 'NS+CF','../figs/NS+CF_Species2.pdf') plotSpecies2(curvesNS_CD, legendtup, 'NS+CD','../figs/NS+CD_Species2.pdf') + +# In[ ]: + + + +