!date
Sat Apr 6 22:36:48 JST 2019
import numpy as np
import pandas as pd
from pylab import meshgrid
from scipy import exp,optimize,log,floor
from scipy.signal import argrelextrema
from scipy.optimize import fsolve
from scipy.integrate import ode
backend = 'dopri5'
# Timer
import time
# To store data
import pickle
# Plotting
%matplotlib inline
# Make inline plots raster graphics
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('pdf', 'svg')
# import modules for plotting and data analysis
import matplotlib.pyplot as plt
from matplotlib import gridspec,rc,colors
import matplotlib.ticker as plticker
# Parameters for seaborn plots
import seaborn as sns
clrs = sns.color_palette("Spectral", 6)
def set_plot_style(usetex=False):
sns.set_style('white', {'axes.linewidth': 0.5})
sns.set(style='white', font_scale=1.1,#context='paper',
rc={'xtick.major.size': 6, 'ytick.major.size': 6, 'legend.fontsize': 14,
'text.usetex': usetex, 'font.family': 'serif', 'font.serif': ['Verdana'],
'text.latex.preamble': r"\usepackage{type1cm}"})
plt.rcParams['xtick.major.size'] = 6
plt.rcParams['xtick.major.width'] = 1
plt.rcParams['ytick.major.size'] = 6
plt.rcParams['ytick.major.width'] = 1
plt.rcParams['xtick.bottom'] = True
plt.rcParams['ytick.left'] = True
set_plot_style(True)
# some shapes
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
clrs = sns.color_palette("colorblind",8)
lwd = 2 # line width by default
# Heaviside function #
Heaviside = lambda x: 1 * (x >= 0)
# Cut function #
cut_func = lambda x: 1. if x > 1. else x
import warnings
warnings.filterwarnings("ignore")
import sys
print(sys.version)
from IPython.core.display import Image, display
!mkdir -p ../figures
!mkdir -p ../figures/pngs
3.7.3 (default, Mar 27 2019, 16:54:48) [Clang 4.0.1 (tags/RELEASE_401/final)]
We have the following dynamics to describe the relative frequency $x$ of drug-tolerant cells:
$$\frac{\mathrm dn(t)}{\mathrm dt} = b\left(1-\frac{\chi(1-x(t))}{e^{\kappa\Delta E}+1}-cx(t)\right)n(t)-dn(t)\,,$$$$\frac{\mathrm dx(t)}{\mathrm dt} = b\left(\frac\chi{e^{\kappa\Delta E}+1}-c\right)x(t)(1-x(t))+\frac{\mu(1-x(t))}{e^{\kappa\Delta E}+1}-\bar\mu e^{-\kappa E^{-}}x(t)\,.$$Dynamics incorporated in Python:
α = 0.3; θ = 0.45
κ = 40.0 # Robustness parameter for the main pathway
Reduction and translocation factors to the expression of the main pathway due to the treatment
A = lambda σ: 1-σ*(1-θ)+1e-6
Production function as a step-like function and corresponding potential function
f = lambda y, σ: A(σ)*(α+(1-α)*Heaviside(y-θ))
U = lambda y, σ: -A(σ)*(α+(1-α)*Heaviside(y-θ))*(y-θ)+(y**2-θ**2)/2
Corresponding potential bariers (will be required for the dynamics)
Eplus = lambda σ: U(θ,σ)-U(f(A(σ),σ),σ)
Eminus = lambda σ: U(θ,σ)-U(f(α*A(σ),σ),σ)
Difference in potential bariers
ΔE = lambda σ: Eplus(σ)-Eminus(σ)
d = 0.13 # death rate per day
b = (0.1*(exp(κ*ΔE(1))+1)-0.14*(exp(κ*ΔE(0))+1))/(exp(κ*ΔE(1))-exp(κ*ΔE(0)))
χ = 1-(0.14*(exp(κ*ΔE(0))+1)-b*exp(κ*ΔE(0)))/b
print("Birth rate: %.4f" % b)
print("Penalty χ: %.4f" % χ)
Birth rate: 0.1402 Penalty χ: 0.3260
ε = 0.01
c_relative = 0.1
c = c_relative*(b-d)/b+(1-c_relative)*χ/(exp(κ*ΔE(0))+1) # cost of resistance
print("Cost of resistance: %.4f" % (c))
print("Cost of resistance, multiplied on b: %.4f" % (b*c))
print("Value of penalized fitness (b*(1-c)-d): %.4f" % (b*(1-c)-d))
Cost of resistance: 0.0083 Cost of resistance, multiplied on b: 0.0012 Value of penalized fitness (b*(1-c)-d): 0.0090
# Transition rates
μ = 1/10.5; μbar = 1/14.
(μ,μbar)
(0.09523809523809523, 0.07142857142857142)
# Time horizon
tmax = 6*30
All parts of Figure 1 were assembled in external program Adobe Illustrator (licence of Hokkaido University).
!convert -flatten -density 150 -trim ../figures/Fig1.pdf -quality 100 ../figures/pngs/Fig1.png
display(Image("../figures/pngs/Fig1.png", width="75%"))
# X = (X[1],X[2]) = (x,n)
ode_rhs = lambda t, X, σ: [b*(χ/(exp(κ*ΔE(σ))+1)-c)*X[0]*(1-X[0])+μ*(1-X[0])/(exp(κ*ΔE(σ))+1)-μbar*exp(-κ*Eminus(σ))*X[0],\
b*(1-χ*(1-X[0])/(exp(κ*ΔE(σ))+1)-c*X[0])*X[1]-d*X[1]]
# Only to track the dynamics of x (relative fraction of resistant cells)
ode_rhs_only_x = lambda t, x, σ: b*(χ/(exp(κ*ΔE(σ))+1)-c)*x*(1-x)+μ*(1-x)/(exp(κ*ΔE(σ))+1)-μbar*exp(-κ*Eminus(σ))*x
plt.rcParams['figure.figsize'] = (4.75, 4.75)
fig, ax = plt.subplots()
ax.margins(0.0)
σs = np.arange(0,1.2,.2)
clrs = sns.color_palette("colorblind",len(σs))
plt.plot([0,1],[0,1],lw=.8,ls="--",color="gray")
for i in range(len(σs)):
σ = σs[i]
plt.plot([0,θ,θ,1],[f(0,σ),f(0,σ),f(1,σ),f(1,σ)],color=clrs[i],lw=lwd,label="%.1f"%σ);
plt.xlabel('Pathway activity level',labelpad=20);
plt.ylabel('Production function',labelpad=12);
legend = plt.legend(loc=2,ncol=1,frameon=1,title='Treatm. intensity',fontsize=11);
plt.setp(legend.get_title(),fontsize=12)
legend.draw_frame(False)
plt.tight_layout()
plt.savefig("../figures/draft/FigS1-A.pdf",format='pdf',bbox_inches='tight')
plt.rcParams['figure.figsize'] = (4.5, 4.5)
fig, ax = plt.subplots()
ax.margins(0.0)
plt.plot([0,1],[0,0],lw=.8,ls="--",color="gray")
for i in range(len(σs)):
σ = σs[i]
yi = np.concatenate([np.linspace(0.,θ,40),np.linspace(θ,1.,40)])
plt.plot(yi,U(yi,σ),color=clrs[i],lw=2,label=σ);
plt.xlabel('Pathway activity level',labelpad=20);
plt.ylabel('Potential function',labelpad=12);
plt.savefig("../figures/draft/FigS1-B.pdf",format='pdf',bbox_inches='tight')
%%capture
%%bash
pdflatex <<TeXScript
\documentclass{standalone}
\usepackage{tikz,graphicx}
\usepackage[T1]{fontenc}
\begin{document}
\begin{tikzpicture}
\node[anchor=north east] {\includegraphics[width=2.95in]{../figures/draft/FigS1-A.pdf}};
\node[anchor=north west] at (.1in,0in) {\includegraphics[width=3.1in]{../figures/draft/FigS1-B.pdf}};
\node[anchor=south east] at (-2.85in,-.04in) {\large {\bf A}};
\node[anchor=south west] at (.06in,-.04in) {\large {\bf B}};
\end{tikzpicture}
\end{document}
TeXScript
mv texput.pdf ../figures/FigS1.pdf
rm texput.*
!convert -flatten -density 150 -trim ../figures/FigS1.pdf -quality 100 ../figures/pngs/FigS1.png
display(Image("../figures/pngs/FigS1.png", width="75%"))
Let's look how $\Delta E$, $\frac{\chi}{e^{\kappa \Delta E}+1}$ and $e^{-\kappa E^{-}}$ depend on $\sigma$
plt.rcParams['figure.figsize'] = (9.5, 3.8)
ax = plt.subplot(121)
plt.subplots_adjust(wspace = 0.3)
plt.plot([0,1],[0,0],lw=.8,ls="--",color="gray")
σi = np.linspace(0.,1.,50);
plt.plot(σi,ΔE(σi),color="k",lw=2,label=r'$\Delta E$');
plt.xlabel('Treatment intensity',labelpad=20);
plt.legend(loc=1,frameon=0);
plt.margins(0.0)
ax.yaxis.set_major_locator(plticker.MultipleLocator(.05))
ax = plt.subplot(122)
plt.plot(σi,χ/(exp(κ*ΔE(σi))+1.),color="k",lw=2,label=r'${\chi}/{(e^{\kappa\Delta E}+1)}$');
plt.plot(σi,exp(-κ*Eminus(σi)),color="k",lw=2,ls='--',label=r'$e^{-\kappa E^{-}}$');
plt.plot(σi,[c_relative]*len(σi),color="k",lw=2,ls='-.',label=r'$c_{\mathrm{relative}}$');
plt.xlabel('Treatment intensity',labelpad=20);
legend = plt.legend(loc=1,frameon=0,ncol=1)
plt.margins(0.0)
ax.yaxis.set_major_locator(plticker.MultipleLocator(.1))
plt.savefig("../figures/FigS2.pdf",format='pdf',bbox_inches='tight')
plt.rcParams['figure.figsize'] = (3.75, 4)
def func_threshold(τ,ε,σ,threshold):
solver = ode(ode_rhs).set_integrator(backend)
solver.set_initial_value([ε,1.]).set_f_params(σ)
solver.integrate(τ)
return solver.y[1]-threshold
σi = [1.]; δσ = .01
t_sol = 450.; t_sol10 = 600.; t_thr = []; t_thr10 = []
while σi[-1]>=0:
t_sol = fsolve(func_threshold,t_sol,args=(ε,σi[-1],1.))[0]; t_thr += [t_sol]
t_sol10 = fsolve(func_threshold,t_sol10,args=(ε,σi[-1],2))[0]; t_thr10 += [t_sol10]
σi += [σi[-1]-δσ]
σi = σi[:-1]
plt.plot(σi,[x+1 for x in t_thr],color="k",lw=lwd,ls='-')
plt.plot(σi,[x+1 for x in t_thr10],color="k",lw=lwd,ls='--')
imx = np.argmax(t_thr)
print("Best outcome: %.1f months at σ = %.2f" % ((1+t_thr[imx])/30,σi[imx]))
print("MTD gives the relapse in %.1f months" % (t_thr[0]/30))
plt.xlabel('Treatment intensity',labelpad=20);
plt.ylabel('Time (months)',labelpad=12)
plt.xlim([0,1])
plt.ylim([0,6*30])
plt.xticks(np.arange(0,1.2,.2))
plt.yticks(np.arange(0,12*30,2*30),np.arange(0,12,2))
plt.margins(0.0)
plt.savefig("../figures/draft/FigS3-A.pdf",format='pdf',bbox_inches='tight')
Best outcome: 2.3 months at σ = 0.64 MTD gives the relapse in 1.8 months
def func_sol(τ,ε,σ):
solver = ode(ode_rhs).set_integrator(backend)
solver.set_initial_value([ε,1.]).set_f_params(σ)
solver.integrate(τ)
return solver.y[1]
σi = [1.]; δσ = .025
tmins = []; reduction = []
while σi[-1]>=0:
τs = np.linspace(0,80,200+1)
sol = [func_sol(τ,ε,σi[-1]) for τ in τs]
index_min = np.argmin(sol)
tmins += [τs[index_min]]
reduction += [1-sol[index_min]]
σi += [σi[-1]-δσ]
σi = σi[:-1]
plt.rcParams['figure.figsize'] = (3.75, 3.25)
plt.plot(σi,[100*x for x in reduction],color="k",lw=lwd,ls='-')
print("Reduction for MTD (sigma = %.2f): %.1f%%"%(σi[0],100*reduction[0]))
plt.xlabel('Treatment intensity',labelpad=20);
plt.ylabel(r'Maximal tumour shrinkage (\%)',labelpad=12)
plt.xticks(np.arange(0,1.2,.2)), plt.yticks(np.arange(0,60,10))
plt.xlim([0,1]), plt.ylim([0,20])
plt.margins(0.0)
plt.savefig("../figures/draft/FigS3-B.pdf",format='pdf',bbox_inches='tight')
Reduction for MTD (sigma = 1.00): 16.8%
plt.rcParams['figure.figsize'] = (4.5, 3.5)
def func_fifty(τ,ε,σ):
solver = ode(ode_rhs_only_x).set_integrator(backend)
solver.set_initial_value([ε]).set_f_params(σ)
solver.integrate(τ)
return solver.y[0]-.5
ε0 = 0.01
σi = [1.];
δσ = .01
t_sol = 50.; t_fifty = [];
while t_sol <= tmax:
t_sol = fsolve(func_fifty,t_sol,args=(ε0,σi[-1]))[0]
t_fifty += [t_sol]
σi += [σi[-1]-δσ]
σi = σi[:-1]
plt.plot(σi,t_fifty,color="k",lw=lwd,ls='-')
# in reverse direction: from 99% resistance toward 50%
σi = [.0];
t_sol = 50.; t_fifty = [];
while t_sol <= tmax:
t_sol = fsolve(func_fifty,t_sol,args=(1-ε0,σi[-1]))[0]
t_fifty += [t_sol]
σi += [σi[-1]+δσ]
σi = σi[:-1]
plt.plot(σi,t_fifty,color="k",lw=lwd,ls='--')
plt.xlabel('Treatment intensity',labelpad=20);
plt.ylabel('Time (days)',labelpad=12)
plt.xlim([0,1]), plt.ylim([0,80])
plt.xticks(np.arange(0,1.2,.2))#, plt.yticks(np.arange(0,tmax/2+60,60))
plt.margins(0.0)
plt.savefig("../figures/draft/FigS3-C.pdf",format='pdf',bbox_inches='tight')
(Caption) Time necessary to acquire 50%-level of resistance starting from initial 0% (solid) vs time necessary to lower the resistance level from 100% to 50% (dashed).
%%capture
%%bash
pdflatex <<TeXScript
\documentclass{standalone}
\usepackage{tikz,graphicx}
\usepackage[T1]{fontenc}
\begin{document}
\begin{tikzpicture}
\node[anchor=north east] {\includegraphics[width=2.95in]{../figures/draft/FigS3-A.pdf}};
\node[anchor=north west] at (.1in,-.125in) {\includegraphics[width=3.1in]{../figures/draft/FigS3-B.pdf}};
\node[anchor=north west] at (-1.6in,-3.3in) {\includegraphics[width=3.1in]{../figures/draft/FigS3-C.pdf}};
\node[anchor=south east] at (-2.85in,-.04in) {\large {\bf A}};
\node[anchor=south west] at (.06in,-.04in) {\large {\bf B}};
\node[anchor=south west] at (-1.65in,-3.4in) {\large {\bf C}};
\end{tikzpicture}
\end{document}
TeXScript
mv texput.pdf ../figures/FigS3.pdf
rm texput.*
# Final result
!convert -flatten -density 150 -trim ../figures/FigS3.pdf -quality 100 ../figures/pngs/FigS3.png
display(Image("../figures/pngs/FigS3.png", width="75%"))
convert: profile 'icc': 'RGB ': RGB color space not permitted on grayscale PNG `../figures/pngs/FigS3.png' @ warning/png.c/MagickPNGWarningHandler/1747.
ymx = 3.4
ymn = .6
ti = np.linspace(0,tmax,1e2+1)
plt.rcParams['figure.figsize'] = (10, 4)
gs = gridspec.GridSpec(1,2,width_ratios=[1,1])
plt.subplot(gs[0])
plt.subplots_adjust(wspace = 0.2)
print("Static treatment")
_i = 0
for σ in σs:
solver = ode(ode_rhs).set_integrator(backend)
solver.set_initial_value([ε,1.]).set_f_params(.5*σ)
sol = []; _k = 0;
while solver.t < max(ti):
solver.integrate(ti[_k])
sol.append([solver.t]+list(solver.y))
_k += 1
plt.plot([x[0] for x in sol],[x[2] for x in sol],color=clrs[_i],lw=2,label="%.1f"%(.5*σ));
print("σ = %.1f leads to %.2f-fold change after %d days (final resistance level = %.1f%%)" %
(.5*σ,sol[-1][2],sol[-1][0],sol[-1][1]*1e2))
_i += 1;
plt.plot([0,max(ti)],[1]*2,lw=.5,ls="-",color="gray",zorder=-2)
plt.ylabel('Fold change in tumor size',labelpad=12);
plt.xlim([0,tmax])
plt.ylim([ymn,ymx])
xtks = np.arange(0,tmax+30,30)
plt.xticks(xtks,[int(x/30) for x in xtks])
plt.xlabel('Time (months)',labelpad=20)
legend = plt.legend(ncol=2,frameon=0,title=r'Applied treatment $\sigma$',loc=2,fontsize=11);
plt.setp(legend.get_title(),fontsize=12)
plt.title('Static treatment', y=1.05)
plt.margins(0.0)
plt.subplot(gs[1])
print("Periodic treatment")
treatment_periodicity = 14 # in days
_i = 0
for σ in σs:
solver = ode(ode_rhs).set_integrator(backend)
solver.set_initial_value([ε,1.])
sol = []; _k = 0;
tme = 0; treatment = True
while (solver.t < max(ti)):
while (tme+treatment_periodicity < ti[_k]):
tme += treatment_periodicity
solver.set_f_params(int(treatment)*σ).integrate(tme)
treatment = not treatment
solver.set_f_params(int(treatment)*σ).integrate(ti[_k])
sol.append([solver.t]+list(solver.y))
_k += 1
plt.plot([x[0] for x in sol],[x[2] for x in sol],color=clrs[_i],lw=2,label="%.1f"%(.5*σ));
print("σ = %.1f leads to %.2f-fold change after %d days (final resistance level = %.1f%%)" %
(σ,sol[-1][2],sol[-1][0],sol[-1][1]*1e2))
_i += 1;
plt.plot([0,max(ti)],[1]*2,lw=.5,ls="-",color="grey",zorder=-2)
plt.xlim([0,tmax])
plt.ylim([ymn,ymx])
xtks = np.arange(0,tmax+30,30)
plt.xticks(xtks,[int(x/30) for x in xtks])
plt.xlabel('Time (months)',labelpad=20)
plt.title('Periodic treatment', y=1.05)
plt.margins(0.0)
legend = plt.legend(ncol=2,frameon=0,title=r'Average treatm. intensity $\bar\sigma$',loc=2,fontsize=11);
plt.setp(legend.get_title(),fontsize=12)
plt.savefig("../figures/draft/Fig4-AB.pdf",format='pdf',bbox_inches='tight')
Static treatment σ = 0.0 leads to 6.04-fold change after 180 days (final resistance level = 0.7%) σ = 0.1 leads to 5.60-fold change after 180 days (final resistance level = 2.8%) σ = 0.2 leads to 4.56-fold change after 180 days (final resistance level = 9.4%) σ = 0.3 leads to 3.07-fold change after 180 days (final resistance level = 24.8%) σ = 0.4 leads to 1.97-fold change after 180 days (final resistance level = 47.7%) σ = 0.5 leads to 1.55-fold change after 180 days (final resistance level = 66.8%) Periodic treatment σ = 0.0 leads to 6.04-fold change after 180 days (final resistance level = 0.7%) σ = 0.2 leads to 5.17-fold change after 180 days (final resistance level = 5.7%) σ = 0.4 leads to 2.75-fold change after 180 days (final resistance level = 30.1%) σ = 0.6 leads to 1.52-fold change after 180 days (final resistance level = 61.4%) σ = 0.8 leads to 1.32-fold change after 180 days (final resistance level = 75.0%) σ = 1.0 leads to 1.34-fold change after 180 days (final resistance level = 80.4%)
plt.rcParams['figure.figsize'] = (5, 3.5)
clr0 = sns.color_palette("Dark2",3)[1]
clrs1 = sns.color_palette("Set1",3)
clrs2 = sns.color_palette("Set1",2)
T = tmax # six months
ti = np.linspace(0,T,1e3+1)
σi = np.linspace(0.,1.,101)
szs = []
treatment_periodicity = 7
for σ in σi:
solver = ode(ode_rhs).set_integrator(backend)
solver.set_initial_value([ε,1.],0.)
tme = 0; treatment = True
while (tme+treatment_periodicity < T):
tme += treatment_periodicity
solver.set_f_params(int(treatment)*σ).integrate(tme)
treatment = not treatment
solver.set_f_params(int(treatment)*σ).integrate(T)
szs = szs+[solver.y[1]]
plt.plot(.5*σi,szs,lw=lwd,ls="-",color=clrs1[1],label='periodic',zorder=2);
res = min(enumerate(szs), key=lambda p: p[1])
print("Best outcome for periodic treatment: sigma = %.2f, fold change = %.4f"%(.5*σi[res[0]],szs[res[0]]))
szs = []
for σ in σi:
solver = ode(ode_rhs).set_integrator(backend)
solver.set_initial_value([ε,1.],0.).set_f_params(1.*σ).integrate(T)
szs = szs+[solver.y[1]]
plt.plot(σi,szs,lw=lwd,ls="-",color=clrs1[0],label='static',zorder=-1);
xp = 10
print("Outcome for static particular treatment: sigma = %.2f, fold change = %.4f"%(σi[xp],szs[xp]))
print("Outcome for static MTD treatment: sigma = %.2f, fold change = %.4f"%(σi[len(σi)-1],szs[len(σi)-1]))
res = min(enumerate(szs), key=lambda p: p[1])
print("Best outcome for static treatment: sigma = %.2f, fold change = %.4f"%(σi[res[0]],szs[res[0]]))
plt.xlabel(r'Average treatment intensity $\bar\sigma$',labelpad=20);
plt.ylabel('Fold change in tumor size',labelpad=12);
plt.xlim([0,1]);
plt.xticks(np.arange(0,1.25,.25));
plt.ylim([0.8,4.2]);
plt.margins(0.0);
legend = plt.legend(frameon=0,title='Regimen',loc=1,fontsize=12); #,fontsize='small'
plt.setp(legend.get_title(),fontsize=12);
# This was obtained in solution of the optimal control problem #
plt.plot([0,1],[1]*2,lw=.8,color=clr0,ls="--",zorder=-4,dashes=(10,5));
plt.savefig("../figures/draft/Fig4-C.pdf",format='pdf',bbox_inches='tight')
Best outcome for periodic treatment: sigma = 0.43, fold change = 1.3379 Outcome for static particular treatment: sigma = 0.10, fold change = 5.5977 Outcome for static MTD treatment: sigma = 1.00, fold change = 2.1330 Best outcome for static treatment: sigma = 0.57, fold change = 1.4832
%%capture
%%bash
pdflatex <<TeXScript
\documentclass{standalone}
\usepackage{tikz,graphicx}
\usepackage[T1]{fontenc}
\begin{document}
\begin{tikzpicture}
\node[anchor=north west] {\includegraphics[width=5.8in]{../figures/draft/Fig4-AB.pdf}};
\node[anchor=north east] at (4.7in,-3.1in) {\includegraphics[width=3.1in]{../figures/draft/Fig4-C.pdf}};
\node[anchor=south east] at (.1in,-.08in) {\large {\bf A}};
\node[anchor=south west] at (3in,-.08in) {\large {\bf B}};
\node[anchor=south west] at (1.45in,-3.18in) {\large {\bf C}};
\end{tikzpicture}
\end{document}
TeXScript
mv texput.pdf ../figures/Fig4.pdf
rm texput.*
# Final result
!convert -flatten -density 150 -trim ../figures/Fig4.pdf -quality 100 ../figures/pngs/Fig4.png
display(Image("../figures/pngs/Fig4.png", width="75%"))
i = 0
prds = [30,36]
for σFixed in [.4,.8]:
plt.rcParams['figure.figsize'] = (3.6, 6.5)
plt.subplot(211)
plt.subplots_adjust(hspace = 0.45)
periods = np.arange(1,T/2,.05)
szs = []; sgms_x = []; resistance_x = [];
for treatment_periodicity in periods:
solver = ode(ode_rhs).set_integrator(backend).set_initial_value([ε,1.])
tme = 0; treatment = True
while (tme+treatment_periodicity < T):
tme += treatment_periodicity
solver.set_f_params(int(treatment)*σFixed).integrate(tme)
treatment = not treatment
solver.set_f_params(int(treatment)*σFixed).integrate(T)
szs = szs+[solver.y[1]]
solver = ode(ode_rhs).set_integrator(backend)
solver.set_initial_value([ε,1.]).set_f_params(.5*σFixed).integrate(T)
sz_static = solver.y[1]
plt.plot([0,max(periods)],[sz_static]*2,color=clr0,lw=1,ls='--',zorder=-1)
for k in range(2):
treatment_periodicity = prds[k]
solver = ode(ode_rhs).set_integrator(backend).set_initial_value([ε,1.])
tme = 0; treatment = True
while (tme+treatment_periodicity < T):
tme += treatment_periodicity
solver.set_f_params(int(treatment)*σFixed).integrate(tme)
treatment = not treatment
solver.set_f_params(int(treatment)*σFixed).integrate(T)
plt.scatter([treatment_periodicity],[solver.y[1]],facecolor=clrs2[k],lw=.5,zorder=3,s=32)
plt.plot(periods,szs,lw=lwd,ls="-",color="k");
plt.title('Applied treatment intensity: %.1f' % (σFixed), y=1.08)
plt.xlabel('Length of drug holidays (days)');
plt.ylabel('Fold change in tumor size',labelpad=12);
plt.xlim([0,max(periods)]); plt.ylim([1,3] if i==1 else [2,5]);
plt.xticks(np.arange(0,T/2+30,30))
plt.margins(0.0)
plt.subplot(212)
szs = []; sgms_x = []; resistance_x = [];
for treatment_periodicity in periods:
sgms_x = sgms_x+[σFixed]
solver = ode(ode_rhs).set_integrator(backend).set_initial_value([ε,1.])
tme = 0; treatment = True
while (tme+treatment_periodicity < T):
tme += treatment_periodicity
solver.set_f_params(int(treatment)*σFixed).integrate(tme)
treatment = not treatment
solver.set_f_params(int(treatment)*σFixed).integrate(T)
szs = szs+[solver.y[1]]
resistance_x = resistance_x+[solver.y[0]]
solver = ode(ode_rhs).set_integrator(backend)
solver.set_initial_value([ε,1.]).set_f_params(.5*σFixed).integrate(T)
sz_constant = solver.y[1]
resistance_constant = solver.y[0]
plt.plot([0,max(periods)],[resistance_constant]*2,color=clr0,lw=1,zorder=-1,ls="--")
for k in range(2):
treatment_periodicity = prds[k]
solver = ode(ode_rhs).set_integrator(backend).set_initial_value([ε,1.])
tme = 0; treatment = True
while (tme+treatment_periodicity < T):
tme += treatment_periodicity
solver.set_f_params(int(treatment)*σFixed).integrate(tme)
treatment = not treatment
solver.set_f_params(int(treatment)*σFixed).integrate(T)
plt.scatter([treatment_periodicity],[solver.y[0]],facecolor=clrs2[k],lw=.5,zorder=3,s=32)
plt.plot(periods,resistance_x,lw=lwd,ls="-",color="k");
plt.xlabel('Length of drug holidays (days)');
plt.ylabel('Final resistance level',labelpad=12);
plt.xlim([0,max(periods)]); plt.ylim([0,1]);
plt.xticks(np.arange(0,T/2+30,30))
plt.margins(0.0)
plt.savefig("../figures/draft/Fig5-CD%d.pdf" % (i+1),format='pdf',bbox_inches='tight')
plt.show()
i += 1
plt.rcParams['figure.figsize'] = (5.2, 1.1)
plt.subplots_adjust(hspace = 0.3)
k = 0; σFixed = .8; n = 3
for treatment_periodicity in [T/(2*n),T/(2*n-1)]:
patches = []
σx = σFixed*((2*n-1)/(2*n) if k==1 else 1)
for i in range(n+1):
plt.plot([2*i*treatment_periodicity,(2*i+1)*treatment_periodicity],[σx]*2,lw=lwd,color=clrs2[k],zorder=-3+2*k)
for l in range(2):
plt.plot([(2*i+l)*treatment_periodicity,(2*i+l)*treatment_periodicity],[σx,0],lw=.75,ls='-',
color=clrs2[k],zorder=-3+2*k)
polygon = Polygon([(2*i*treatment_periodicity,0),(2*i*treatment_periodicity,σx),
((2*i+1)*treatment_periodicity,σx),((2*i+1)*treatment_periodicity,0)], True)
patches.append(polygon)
p = PatchCollection(patches,alpha=.5,color=clrs2[k], zorder=-4+2*k);
plt.gca().add_collection(p)
k += 1
plt.xlim([0,T]); plt.ylim([0,.9])
xtks = np.arange(0,T+30,30); plt.xticks(xtks,[x if k else "" for x in xtks])
plt.ylabel('Treatment',labelpad=16);
plt.xlabel('Time (days)',labelpad=10);
ytks = [0,σx,σFixed]; plt.yticks([]);
sns.despine(offset=0)
plt.margins(0.0)
plt.savefig("../figures/draft/Fig5-B.pdf",format='pdf',bbox_inches='tight')
# Change to True if you need recalculation
recalc3 = False
def func_σ(treatment_periodicity,σx,out):
solver = ode(ode_rhs).set_integrator(backend).set_initial_value([ε,1.])
tme = 0; treatment = True
if out: print(σx)
while (tme+treatment_periodicity < T):
tme += treatment_periodicity
solver.set_f_params(int(treatment)*σx).integrate(tme)
treatment = not treatment
solver.set_f_params(int(treatment)*σx).integrate(T)
return solver.y
Difference:
step_periods = 1; step_σFixed = .01
Periods, σFixed = np.meshgrid(np.arange(1,T/2+step_periods,step_periods),np.arange(0,1.+step_σFixed,step_σFixed))
periods_n, σFixed_n = Periods.shape
from pathlib import Path
datafile = '../figures/draft/FigS4.pkl'
if (recalc3|(not Path(datafile).is_file())):
time0=time.time()
Resistance = Periods*0
FoldChange = Periods*0
for xk in range(periods_n):
for yk in range(σFixed_n):
Resistance[xk,yk], FoldChange[xk,yk] = func_σ(Periods[xk,yk],σFixed[xk,yk],False)
print("This proccess took %0.1f minutes" % ((time.time()-time0)/60.))
# Saving the data
file_to_write = open(datafile,'wb')
pickle.dump([Resistance, FoldChange], file_to_write)
file_to_write.close()
recalc3 = False
else:
# Loading the data
file_to_load = open(datafile,'rb')
Resistance, FoldChange = pickle.load(file_to_load)
file_to_load.close()
This proccess took 8.2 minutes
xmx = T/2
plt.rcParams['figure.figsize'] = (12, 4)
plt.subplot(121)
plt.subplots_adjust(wspace = 0.125)
lvls = [1.1,1.2,1.3,1.4,1.5,1.6,1.8,2,2.5,3,3.5,4,5,7]
tks = [1.1,1.2,1.4,1.6,2,3,4,5,7]
plt.rcParams['figure.figsize'] = (8, 5)
plt.contourf(Periods,σFixed,FoldChange,cmap="Spectral_r",zorder=-2,alpha=.7,levels=lvls)
cbar = plt.colorbar(pad=.06,aspect=18)
plt.xticks(np.arange(0,xmx+15,15));
sz=20; zmax = 6.
cbar.set_ticks(tks)
pts_x = []; pts_y = []; pts_z = []
for n in range(len(Periods[1,:])):
res = min(enumerate(FoldChange[:,n]), key=lambda p: p[1])
pts_x = np.append(pts_x,Periods[1,n])
pts_y = np.append(pts_y,σFixed[res[0],n])
pts_z = np.append(pts_z,res[1])
res = min(enumerate(pts_z), key=lambda p: p[1])
plt.scatter([pts_x[res[0]]],[pts_y[res[0]]],facecolor="red",marker=(5,1),s=80,linewidth=.5,color="black",zorder=4);
print("Best outcome: period = %.2f, sigma = %.2f, fold change = %.4f"%(pts_x[res[0]],pts_y[res[0]],pts_z[res[0]]))
print("Nearest point: period = %.2f, sigma = %.2f, fold change = %.4f"%(pts_x[res[0]-1],pts_y[res[0]-1],pts_z[res[0]-1]))
plt.plot(pts_x,pts_y,lw=lwd,ls=":",color="black",zorder=3);
plt.xlabel('Length of drug holidays (days)',labelpad=20);
plt.ylabel('Treatment intensity',labelpad=12);
plt.yticks(np.arange(0,1.2,.2))
plt.xlim([0,xmx]); plt.ylim([0,1])
plt.title(r'Fold change', y=1.06)
plt.subplot(122)
plt.contourf(Periods,σFixed,Resistance,cmap="RdYlBu_r",
levels=np.arange(0,1.1,.1),zorder=-2,alpha=.7)
cbar = plt.colorbar(pad=.06,aspect=18)
cbar.set_ticks(np.arange(0,1.2,.2))
cbar.set_ticklabels([int(x*100) for x in np.arange(0,1.2,.2)])
plt.title('Resistance level (\%)', y=1.08)
plt.scatter([pts_x[res[0]]],[pts_y[res[0]]],facecolor="red",marker=(5,1),
s=60,lw=.5,color="black",zorder=4);
plt.plot(pts_x,pts_y,lw=lwd,ls=":",color="black",zorder=3);
plt.xlabel('Length of drug holidays (days)',labelpad=20);
plt.xticks(np.arange(0,xmx+15,15))
plt.yticks(np.arange(0,1.2,.2))
plt.xlim([0,xmx]); plt.ylim([0,1])
plt.savefig("../figures/draft/FigS4-AB.pdf",format='pdf',bbox_inches='tight')
Best outcome: period = 4.00, sigma = 0.88, fold change = 1.2998 Nearest point: period = 3.00, sigma = 0.88, fold change = 1.3316
%%capture
%%bash
pdflatex <<TeXScript
\documentclass{standalone}
\usepackage{tikz,graphicx}
\usepackage[T1]{fontenc}
\begin{document}
\begin{tikzpicture}
\node[anchor=south] {\includegraphics[width=6.2in]{../figures/draft/FigS4-AB.pdf}};
\node[anchor=south east] at (-2.93in,2.66in) {\large {\bf A}};
\node[anchor=south east] at (.25in,2.66in) {\large {\bf B}};
\end{tikzpicture}
\end{document}
TeXScript
mv texput.pdf ../figures/FigS4.pdf
rm texput.*
# Final result
!convert -flatten -density 150 -trim ../figures/FigS4.pdf -quality 100 ../figures/pngs/FigS4.png
display(Image("../figures/pngs/FigS4.png", width="75%"))
** ¡It takes quite a substantial amount of time! **
# Change to True if you need recalculation
recalc2 = False
def funcX(treatment_cycle,x,σ_x):
solver = ode(ode_rhs).set_integrator(backend).set_initial_value([ε,1.])
tme = 0
while (tme+cycle < T):
tme += treatment_cycle*(1+x)/2
solver.set_f_params(σ_x).integrate(tme)
tme += treatment_cycle*(1-x)/2
solver.set_f_params(0.).integrate(tme)
if (tme+cycle*(1+x)/2 < T):
tme += treatment_cycle*(1+x)/2
solver.set_f_params(σ_x).integrate(tme)
solver.set_f_params(0.).integrate(T)
else:
solver.set_f_params(σ_x).integrate(T)
return solver.y
##########
cycle = 8
##########
step_x = .01; step_σ = .01
Xs, σs = np.meshgrid(np.arange(-1.,1.+step_x,step_x),np.arange(0,1.+step_σ,step_σ))
x_n, σ_n = Xs.shape
from pathlib import Path
datafile2 = "../figures/draft/FigS5.pkl"
if (recalc2|(not Path(datafile2).is_file())):
time0=time.time()
FoldChangeX = Xs*0
ResistanceX = Xs*0
TreatmentX = Xs*0
for xk in range(x_n):
for yk in range(σ_n):
ResistanceX[xk,yk], FoldChangeX[xk,yk] = funcX(cycle,Xs[xk,yk],σs[xk,yk])
TreatmentX[xk,yk] = σs[xk,yk]
print("This proccess took %0.1f minutes" % ((time.time()-time0)/60.))
# Saving the data
file_to_save = open(datafile2,'wb')
pickle.dump([TreatmentX, ResistanceX, FoldChangeX], file_to_save)
file_to_save.close()
recalc2 = False
else:
# Loading the data
file_to_load = open(datafile2,'rb')
TreatmentX, ResistanceX, FoldChangeX = pickle.load(file_to_load)
file_to_load.close()
This proccess took 18.5 minutes
plt.rcParams['figure.figsize'] = (12, 4.5)
plt.subplot(121)
plt.subplots_adjust(hspace = 0.3, wspace = .15)
xlb = 'Fraction of the cycle with drug administration (\%)'
lvls = [1.1,1.2,1.3,1.4,1.5,1.6,1.8,2,2.5,3,3.5,4,5,7]
tks = [1.1,1.2,1.4,1.6,2,3,4,5,7]
cut_factor = 3
pts_x = []; pts_y = []; pts_z = []
for n in range(1,len(Xs[1,:])-1): #the first and the last elements are removed
res = min(enumerate(FoldChangeX[:,n]), key=lambda p: p[1])
pts_x = np.append(pts_x,Xs[1,n])
pts_y = np.append(pts_y,σs[res[0],n])
pts_z = np.append(pts_z,res[1])
res = min(enumerate(pts_z), key=lambda p: p[1])
plt.scatter([pts_x[res[0]]],[pts_y[res[0]]],facecolor="red",marker=(5,1),s=80,linewidth=.5,color="black",zorder=4);
print("Best outcome: ψ = %.2f, σ = %.2f, fold change = %.3f" % ((1+pts_x[res[0]])/2,pts_y[res[0]],pts_z[res[0]]))
plt.contourf(Xs,σs,FoldChangeX,cmap="Spectral_r",zorder=-2,alpha=.7,
levels=lvls)
cbar = plt.colorbar(pad=.06,aspect=18)
xtks = np.arange(-1,1.5,.5)
plt.xticks(xtks,["%d" % int(100*(x+1)/2) for x in xtks])
plt.xlim([-1,1]); plt.ylim([0,1])
cbar.set_ticks(tks)
plt.title('Fold change', y=1.04)
plt.ylabel('Treatment intensity',labelpad=12);
plt.yticks(np.arange(0,1.2,.2))
plt.xlabel(xlb,labelpad=12);
plt.subplot(122)
plt.scatter([pts_x[res[0]]],[pts_y[res[0]]],facecolor="red",marker=(5,1),
s=60,lw=.5,color="black",zorder=4);
# plt.plot(pts_x,pts_y,lw=lwd,ls=":",color="black",zorder=3);
plt.contourf(Xs,σs,ResistanceX,cmap="RdYlBu_r",levels=np.arange(0,1.1,.1),zorder=-2,alpha=.7)
plt.title('Resistance level (\%)', y=1.04)
cbar = plt.colorbar(pad=.06,aspect=18)
cbar.set_ticks(np.arange(0,1.2,.2))
cbar.set_ticklabels([int(x*100) for x in np.arange(0,1.2,.2)])
plt.xlabel(xlb,labelpad=12);
plt.ylabel('Treatment intensity',labelpad=12);
plt.xticks(xtks,["%d" % int(100*(x+1)/2) for x in xtks])
plt.yticks(np.arange(0,1.2,.2))
plt.xlim([-1,1]); plt.ylim([0,1])
plt.savefig("../figures/draft/FigS5-AB.pdf",format='pdf',bbox_inches='tight')
Best outcome: ψ = 0.56, σ = 0.82, fold change = 1.292
plt.rcParams['figure.figsize'] = (4.2, 1.2)
x0 = -0.2
plt.plot([0,(1+x0)/2],[1,1],lw=lwd,color="k")
plt.plot([(1+x0)/2,1],[0,0],lw=lwd,color="k")
plt.plot([(1+x0)/2,(1+x0)/2],[1,0],lw=.75,color="k",ls='-',zorder=-3)
p = PatchCollection([Polygon([(0,0),(0,1),((1+x0)/2,1),((1+x0)/2,0)], True)],alpha=.5,color="grey");
plt.gca().add_collection(p)
xtks = np.arange(0,1.25,.25)
plt.xlim([0,1]); plt.ylim([0,1.2])
plt.xticks(xtks,["%d" % int(100*x) for x in xtks])
plt.title('Particular treatment regime',y=1.06);
plt.xlabel('Fraction of the treatment cycle (\%)',labelpad=8);
#plt.yticks([0,1],[0,r'$\sigma=\min\left(\frac{2\bar\sigma}{1+\psi},1\right)$']);
sns.despine(offset=0)
plt.savefig("../figures/draft/FigS5-C.pdf",format='pdf',bbox_inches='tight')
%%capture
%%bash
pdflatex <<TeXScript
\documentclass{standalone}
\usepackage{tikz,graphicx}
\usepackage[T1]{fontenc}
\begin{document}
\begin{tikzpicture}
\node[anchor=north west] {\includegraphics[width=\textwidth]{../figures/draft/FigS5-AB.pdf}};
\node[anchor=north west] at (.25in,-2.32in) {\includegraphics[width=.36\textwidth]{../figures/draft/FigS5-C.pdf}};
\node[anchor=south east] at (.15in,-.1in) {{\bf A}};
\node[anchor=south east] at (2.6in,-.1in) {{\bf B}};
\node[anchor=south east] at (.15in,-2.42in) {{\bf C}};
\end{tikzpicture}
\end{document}
TeXScript
mv texput.pdf ../figures/FigS5.pdf
rm texput.*
# Final result
!convert -flatten -density 150 -trim ../figures/FigS5.pdf -quality 100 ../figures/pngs/FigS5.png
display(Image("../figures/pngs/FigS5.png", width="75%"))
We generate only output for two treatment and we will use it later in another R-notebook.
def funcXtrj(treatment_cycle,x,σ_x,nsteps=10):
solver = ode(ode_rhs).set_integrator(backend).set_initial_value([ε,1.])
sol = [[0.,σ_x,ε,1.]]; tme = 0
while (tme+treatment_cycle < T):
デルタ = treatment_cycle*(1+x)/2
ti = np.linspace(tme,tme+デルタ,nsteps)
_k = 0
while (solver.t < max(ti)):
solver.set_f_params(σ_x).integrate(ti[_k])
if _k:
sol.append([solver.t,σ_x]+list(solver.y))
_k += 1
tme += デルタ
デルタ = treatment_cycle*(1-x)/2
ti = np.linspace(tme,tme+デルタ,nsteps)
_k = 0
while (solver.t < max(ti)):
solver.set_f_params(0.0).integrate(ti[_k])
if _k:
sol.append([solver.t,0.0]+list(solver.y))
_k += 1
tme += デルタ
if (tme+treatment_cycle*(1+x)/2 < T):
デルタ = treatment_cycle*(1+x)/2
ti = np.linspace(tme,tme+デルタ,nsteps,endpoint=False)
_k = 0
while (solver.t < max(ti)):
solver.set_f_params(σ_x).integrate(ti[_k])
if _k:
sol.append([solver.t,σ_x]+list(solver.y))
_k += 1
tme += デルタ
ti = np.linspace(tme,T,nsteps,endpoint=True)
_k = 0
while (solver.t < max(ti)):
solver.set_f_params(0.0).integrate(ti[_k])
if _k:
sol.append([solver.t,0.0]+list(solver.y))
_k += 1
else:
ti = np.linspace(tme,T,nsteps,endpoint=True)
_k = 0
while (solver.t < max(ti)):
solver.set_f_params(σ_x).integrate(ti[_k])
if _k:
sol.append([solver.t,σ_x]+list(solver.y))
_k += 1
return sol
pd.DataFrame(funcXtrj(8,0,.88,100)).to_csv('../figures/draft/Fig7-trj_periodic.csv',index=False,header=False)
pd.DataFrame(funcXtrj(8,1,.57,100)).to_csv('../figures/draft/Fig7-trj_const.csv',index=False,header=False)
Draft figures 6A, 6B and 6C are obtained in the notebook "C1. Sensitivity analysis.ipynb". Then the figure was assembled in external program Adobe Illustrator (licence of Hokkaido University) and renamed to Fig 8.
# Final result
!convert -flatten -density 150 -trim ../figures/Fig8.pdf -quality 100 ../figures/pngs/Fig8.png
display(Image("../figures/pngs/Fig8.png", width="75%"))
It was obtained in the notebook "D. Field of optimal trajectories.ipynb". Here we modify it a bit using tikz/LaTeX
%%capture
%%bash
pdflatex <<TeXScript
\documentclass{standalone}
\usepackage{tikz,graphicx}
\usepackage[T1]{fontenc}
\begin{document}
\begin{tikzpicture}
\node {\includegraphics[scale=.85]{../figures/draft/Fig2-1.pdf}};
\node at (-1.0in,.86in) {\scalebox{1.2}{$\mathcal S_1$}};
\node at (.36in,.3in) {\scalebox{1.2}{$\mathcal S_2$}};
\node at (-1.4in,.52in) {$\bar x_2$};
\end{tikzpicture}
\end{document}
TeXScript
mv texput.pdf ../figures/Fig3.pdf
rm texput.*
# Final result
!convert -flatten -density 150 -trim ../figures/Fig3.pdf -quality 100 ../figures/pngs/Fig3.png
display(Image("../figures/pngs/Fig3.png", width="45%"))
The constituent parts were obtained in the notebook "D. Field of optimal trajectories.ipynb"
%%capture
%%bash
pdflatex <<TeXScript
\documentclass{standalone}
\usepackage{tikz,graphicx}
\usepackage[T1]{fontenc}
\begin{document}
\begin{tikzpicture}
\node[anchor=south east] at (0,.1in) {\includegraphics[width=2.5in]{../figures/draft/Fig3-A.pdf}};
\node[anchor=south west] at (.25in,.1in) {\includegraphics[width=2.5in]{../figures/draft/Fig3-B.pdf}};
\node[anchor=north east] {\includegraphics[width=2.5in]{../figures/draft/Fig3-C.pdf}};
\node[anchor=north west] at (.25in,0) {\includegraphics[width=2.5in]{../figures/draft/Fig3-D.pdf}};
\node[anchor=south east] at (-2.5in,1.9in) {\large {\bf A}};
\node[anchor=south west] at (.15in,1.9in) {\large {\bf B}};
\node[anchor=south east] at (-2.5in,-.05in) {\large {\bf C}};
\node[anchor=south west] at (.15in,-.05in) {\large {\bf D}};
\node[anchor=south east] at (-.36in,-2.05in) {Treatment intensity $\sigma$};
\node[anchor=south east] at (2.51in,-2.05in) {Treatment intensity $\sigma$};
\end{tikzpicture}
\end{document}
TeXScript
mv texput.pdf ../figures/Fig2.pdf
rm texput.*
# Final result
!convert -flatten -density 150 -trim ../figures/Fig2.pdf -quality 100 ../figures/pngs/Fig2.png
display(Image("../figures/pngs/Fig2.png", width="45%"))
It was obtained in the notebook "E. Figures 6 and 7.ipynb"
It's combined from draft Figure 4 (previous ordering)
%%capture
%%bash
pdflatex <<TeXScript
\documentclass{standalone}
\usepackage{tikz,graphicx}
\usepackage[T1]{fontenc}
\begin{document}
\begin{tikzpicture}
\node[anchor=south] {\includegraphics[scale=1]{../figures/draft/Fig4.pdf}};
\node[anchor=north] at (.36in,.6in) {\(T_0\)};
\node[anchor=north] at (.73in,.6in) {\(T_{-}\)};
\node[anchor=north] at (-1.48in,.6in) {\(T_{+}\)};
\node[anchor=north] at (.76in,1.7in) {(\emph{i})};
\node[anchor=north] at (-1.15in,1.7in) {(\emph{ii})};
\node[anchor=north] at (-1.16in,2.42in) {(\emph{iii})};
\node[anchor=north] at (-.2in,2.3in) {dose-sparing regimen};
\end{tikzpicture}
\end{document}
TeXScript
mv texput.pdf ../figures/Fig6.pdf
rm texput.*
# Final result
!convert -flatten -density 150 -trim ../figures/Fig6.pdf -quality 100 ../figures/pngs/Fig6.png
display(Image("../figures/pngs/Fig6.png", width="45%"))
It was obtained in the notebook "E. Figures 4 and 5.ipynb". Then the figure was assembled in external program Adobe Illustrator (licence of Hokkaido University).
# Final result
!convert -flatten -density 150 -trim ../figures/Fig5.pdf -quality 100 ../figures/pngs/Fig5.png
display(Image("../figures/pngs/Fig5.png", width="45%"))
%%capture
%%bash
pdflatex <<TeXScript
\documentclass{standalone}
\usepackage{tikz,graphicx}
\usepackage[T1]{fontenc}
\begin{document}
\begin{tikzpicture}
\node[anchor=south] {\includegraphics[width=6.2in]{../figures/draft/FigS8.pdf}};
\node[anchor=south east] at (-2.93in,6.2in) {\large {\bf A}};
\node[anchor=south east] at (.35in,6.2in) {\large {\bf B}};
\node[anchor=south east] at (-2.93in,3.1in) {\large {\bf C}};
\node[anchor=south east] at (.35in,3.1in) {\large {\bf D}};
\end{tikzpicture}
\end{document}
TeXScript
mv texput.pdf ../figures/FigS8.pdf
rm texput.*
# Final result
!convert -flatten -density 150 -trim ../figures/FigS8.pdf -quality 100 ../figures/pngs/FigS8.png
display(Image("../figures/pngs/FigS8.png", width="45%"))
%%capture
%%bash
pdflatex <<TeXScript
\documentclass{standalone}
\usepackage{tikz,graphicx}
\usepackage[T1]{fontenc}
\begin{document}
\begin{tikzpicture}
\node[anchor=south east] {\includegraphics{../figures/draft/FigSXa.pdf}};
\node[anchor=south west] {\includegraphics{../figures/draft/FigSXb.pdf}};
\node[anchor=south east] at (-3.5in,3.3in) {\large {\bf A}};
\node[anchor=south east] at (.3in,3.3in) {\large {\bf B}};
\end{tikzpicture}
\end{document}
TeXScript
mv texput.pdf ../figures/FigS6.pdf
rm texput.*
# Final result
!convert -flatten -density 150 -trim ../figures/FigS6.pdf -quality 100 ../figures/pngs/FigS6.png
display(Image("../figures/pngs/FigS6.png", width="45%"))