Page MenuHomec4science
Files F60078979

spent_media_model_interactions.py
No OneTemporary
Actions

File Metadata

Created
Sat, Apr 27, 08:09

spent_media_model_interactions.py
View Options

#!/usr/bin/env python
# coding: utf-8
# 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[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[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_temp<N:
idx_temp += 1
x_temp = x.integrate(x.t+dt)
xout[idx_temp] = x_temp
return t, xout
# In[5]:
def draw_subplot(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="Time [a.u.]",ylabel="Abundance/conc. [a.u.]", title=title_text)
return axtmp
# 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[101]:
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[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"
draw_subplot(ax[1], curves[0], curves[2], legendtup_s2, title_text, 'r',num_C)
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"
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$"
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"
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[237]:
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"
draw_subplot_nolab(ax[0,1], curvesEC[0], curvesEC[2], legendtup_s2, title_text, 'r', num_C)
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"
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$"
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"
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[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[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)]
# In[240]:
#### 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')
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)
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')
AUC
# In[151]:
# 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"]]
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()
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()
fig.savefig('../figs/At_Ct_Oa_AUCs.pdf',format='pdf',dpi=300,bbox_inches='tight')
# In[224]:
# 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
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)
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')
# In[225]:
# 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
fig, axs = plt.subplots(1, len(kpr13_ic), figsize=(1.5 * len(kpr13_ic), 3),
constrained_layout=True, squeeze=False)
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]))
fig.colorbar(psm, ax=axs[0,len(kpr13_ic)-1])
fig.savefig('../figs/sweep_kpr13_comparison.pdf',format='pdf',dpi=300,bbox_inches='tight')
# In[229]:
# 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
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)
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]))
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[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[223]:
# Generate figure to compare curves under different scenarios
curvesNS_CF = all_curves[6]
curvesNS_CD = all_curves[7]
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[ ]:

Event Timeline

c4science · Help