%pylab inline
from scipy.integrate import odeint
import matplotlib.patches as patches
Populating the interactive namespace from numpy and matplotlib
#computes volume time courses of 2 coupled droplets containing the formose reaction
def ODEEXCHP5(VA,VB,vv,PH2O,P1,P2,kk,krel,T,C1A,C2A,C1B,C2B,nplus):
# ODE solver parameters
TT=1
abserr = 1.0e-10
relerr = 1.0e-11
stoptime = T*1.0
numpoints = T*100*nplus
# Create the time samples for the output of the ODE solver.
t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]
k12=kk*krel #C1+C2 aldol rate
k13=kk*krel #C1+C3 aldol rate
k4=kk #C4 cleavage by retro-aldol
k22=kk*krel #C2+C2->C4 aldol
Ninert=0.00
[vvH2O,vvC1,vvC2,vvC3,vvC4]=vv
k=[P1,P2,k12,k13,k4,PH2O,k22,Ninert,vvH2O,vvC1,vvC2,vvC3,vvC4]
NL0=zeros(10)
NL0[0]=VA*(1.0-vvC1*C1A-vvC2*C2A)/vvH2O #VA 90
NL0[1]=VA*C1A #VB 60
NL0[2]=VA*C2A #N2A #3.0 * NL0[0] #N1A
NL0[3]=0.0 #* NL0[0] #N3A
NL0[4]=0.0 #N4A
NL0[5]=VB *(1.0-vvC1*C1B-vvC2*C2B)/vvH2O #NH2OB
NL0[6]=VB * C1B #N1B 3.0*NL0[1] #N1B
NL0[7]=VB * C2B #N2B
NL0[8]=0.0 #N3B
NL0[9]=0.0 #N4B
# Call the ODE solver.
NLF = odeint(vexch7, NL0, t, args=(k,),
atol=abserr, rtol=relerr)
VAF,VBF=[],[]
NNA=[[],[],[],[],[]]
NNB=[[],[],[],[],[]]
for i in range(len(t)):
VAF.append(NLF[int(i)][0]*vvH2O+NLF[int(i)][1]*vvC1+NLF[int(i)][2]*vvC2+NLF[int(i)][3]*vvC3+NLF[int(i)][4]*vvC4)
VBF.append(NLF[int(i)][5]*vvH2O+NLF[int(i)][6]*vvC1+NLF[int(i)][7]*vvC2+NLF[int(i)][8]*vvC3+NLF[int(i)][9]*vvC4)
for k in range(5):
NNA[k].append(NLF[int(i)][k])
NNB[k].append(NLF[int(i)][5+k])
return VAF,VBF,NNA,NNB
#set of equation of the model, computes dX/dt = f(X) where X is the vector of moles in each drop (water and CNs)
def vexch7(NN,t,K):
f=[]
P1,P2=K[0],K[1] #permeabilities of C1 and C2
k12,k13,k4,PH2O,k22=K[2:7] #reaction rates of formose; !!!in the NMR code, Ntot is included in k, not here!!!
Ninert,vH2O,vC1,vC2,vC3,vC4=K[7:13] #molar volume of the species
NH2OA,N1A,N2A,N3A,N4A=NN[0],NN[1],NN[2],NN[3],NN[4] #list of moles in droplet A
NH2OB,N1B,N2B,N3B,N4B=NN[5],NN[6],NN[7],NN[8],NN[9] #list of moles in droplet B
NTA=NH2OA+N1A+N2A+N3A+N4A #total number of moles in droplet A
NTB=NH2OB+N1B+N2B+N3B+N4B #total number of moles in droplet B
VA=NH2OA*vH2O+N1A*vC1+N2A*vC2+N3A*vC3+N4A*vC4 #total volume of droplet A
VB=NH2OB*vH2O+N1B*vC1+N2B*vC2+N3B*vC3+N4B*vC4 #total volume of droplet B
PP=VA**(2.0/3.0)*VB**(2.0/3.0)/(VA**(2.0/3.0)+VB**(2.0/3.0)) #effective exchange area ()
JH2OBA=PH2O*PP*(NH2OA/NTA-NH2OB/NTB) #flux of water in mol/min
JN1BA=P1*PP*(N1A/NTA-N1B/NTB) #flux of C1
JN2BA=P2*PP*(N2A/NTA-N2B/NTB) #flux of C2
#in droplet A
f.append((-JH2OBA)) #H2OA
f.append((-JN1BA - k12*N1A*N2A/NTA - k13*N1A*N3A/NTA)) #variation of C1 due to fluxes and reactions
f.append((-JN2BA - k12*N1A*N2A/NTA + 2.0*k4*N4A -2.0*k22*N2A**2/NTA)) #variation of C2
f.append((-k13*N1A*N3A/NTA + k12*N1A*N2A/NTA)) #variation of C3
f.append((+k13*N1A*N3A/NTA - k4*N4A + k22*N2A**2/NTA)) #variation of C4
#in droplet B
f.append((+JH2OBA)) #H2OB
f.append((+JN1BA - k12*N1B*N2B/NTB - k13*N1B*N3B/NTB))
f.append((+JN2BA - k12*N1B*N2B/NTB + 2.0*k4*N4B -2.0*k22*N2B**2/NTB))
f.append((-k13*N1B*N3B/NTB + k12*N1B*N2B/NTB))
f.append((+k13*N1B*N3B/NTB - k4*N4B + k22*N2B**2/NTB))
return f
Experimentally, a 'flat line' is observed for $C1$ exchange at 8M vs 0.8M. Below, different permeability regimes are illustrated for this scenario, using the above model without reaction. As expected, the toy model gives no discernable growth when $P_{C1}/P_{H_2O} \approx v_{H2O}/v_{C1}$. Flatness is obtained here for $P_{C1}/P_{H_2O}=0.6$. However, when $P_{C1}\gg P_{H_2O}$, the slight initial volume divergence (due to exchange of C1) may occur on a timescale that is not experimentally accessible, so that volume curves would appear flat.
# test different ratios between the permeabilities of C1 and C2
vv=0.018,0.029,0.041,0.053,0.065 #list of molar volumes modified by PN
kk=0.0 #no reaction
krel=0.0 #no reaction
V0=80
Vc=120
P0=30 #done for 20 um for C1, 30=600/20
T=260
C10,C20,C1c,C2c=0.8,0.0,0.08,0 #only C1 exchanges
P1x=0.01*P0
P1a=0.1*P0
P1b=1*P0
P1c=10.0*P0
P2=0.05*P0
[V1q,V2q,_,_]=ODEEXCHP5(V0,Vc,vv,P0,P1x,P2,kk,krel,T,C10,C20,C1c,C2c,1)
[V3q,V4q,_,_]=ODEEXCHP5(V0,Vc,vv,P0,P1a,P2,kk,krel,T,C10,C20,C1c,C2c,1)
[V5q,V6q,_,_]=ODEEXCHP5(V0,Vc,vv,P0,P1b,P2,kk,krel,T,C10,C20,C1c,C2c,1)
[V7q,V8q,_,_]=ODEEXCHP5(V0,Vc,vv,P0,P1c,P2,kk,krel,T,C10,C20,C1c,C2c,1)
P1r=[P1x/P0,P1a/P0,P1b/P0,P1c/P0] #list of P1/P0 ratios used
matplotlib.pyplot.subplots_adjust(left=None, bottom=None, right=1.9, top=2, wspace=0.3, hspace=0.5)
plt.subplot(221)
plt.plot(np.linspace(-20,240,T),V1q[::100],'b',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-20,240,T),V2q[::100],'r',linewidth=3,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,160)
plt.text(20,20,r'$\frac{P_1}{P_{H_2O}}=$'+str(P1r[0]),size=25)
plt.legend(loc=4,fontsize=16)
plt.subplot(222)
plt.plot(np.linspace(-20,240,T),V3q[::100],'b',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-20,240,T),V4q[::100],'r',linewidth=3,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,160)
plt.text(20,20,r'$\frac{P_1}{P_{H_2O}}=$'+str(P1r[1]),size=25)
plt.legend(loc=4,fontsize=16)
plt.subplot(223)
plt.plot(np.linspace(-20,240,T),V5q[::100],'b',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-20,240,T),V6q[::100],'r',linewidth=3,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,160)
plt.text(20,20,r'$\frac{P_1}{P_{H_2O}}=$'+str(P1r[2]),size=25)
plt.legend(loc=4,fontsize=16)
plt.subplot(224)
plt.plot(np.linspace(-20,240,T),V7q[::100],'b',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-20,240,T),V8q[::100],'r',linewidth=3,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,160)
plt.text(20,20,r'$\frac{P_1}{P_{H_2O}}=$'+str(P1r[3]),size=25)
plt.legend(loc=4,fontsize=16)
#plt.savefig('P1overPH20_08vs008.eps',dpi=1000,bbox_inches='tight')
<matplotlib.legend.Legend at 0x27c82d367f0>
Experimentally, we observe that for larger species ($C_n, n \geq 3$), in the absence of the formose reaction, volume time courses appear as linear lines and the intersection point position changes with distance and sugar size (their molar volume). We illustrate these trends with the model below.
We first test the effect of distance for a difference of $C_2$ concentration for different droplet spacings.
vv=0.018,0.029,0.041,0.053,0.065
V0=80
Vc=120
d=[8,14,20,30,40]
P0L=600
T=260
C10,C20,C1c,C2c=0.0,0.8,0.0,0.08
krel=0 #no reaction
kk=0 #no reaction
#Transport coefficients
P1,P2=0.001*P0L,0.001*P0L
[V1q,V2q,_,_]=ODEEXCHP5(V0,Vc,vv,P0L/d[0],P1/d[0],P2/d[0],kk,krel,T,C10,C20,C1c,C2c,1)
[V3q,V4q,_,_]=ODEEXCHP5(V0,Vc,vv,P0L/d[1],P1/d[1],P2/d[1],kk,krel,T,C10,C20,C1c,C2c,1)
[V5q,V6q,_,_]=ODEEXCHP5(V0,Vc,vv,P0L/d[2],P1/d[2],P2/d[2],kk,krel,T,C10,C20,C1c,C2c,1)
[V7q,V8q,_,_]=ODEEXCHP5(V0,Vc,vv,P0L/d[3],P1/d[3],P2/d[3],kk,krel,T,C10,C20,C1c,C2c,1)
[V9q,V10q,_,_]=ODEEXCHP5(V0,Vc,vv,P0L/d[4],P1/d[4],P2/d[4],kk,krel,T,C10,C20,C1c,C2c,1)
lbl=d
matplotlib.pyplot.subplots_adjust(left=None, bottom=None, right=1.9, top=3, wspace=0.3, hspace=0.5)
plt.subplot(321)
plt.plot(np.linspace(-20,240,T),V1q[::100],'b',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-20,240,T),V2q[::100],'r',linewidth=3,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,160)
plt.text(20,20,r'$d=$'+str(lbl[0]),size=25)
plt.legend(loc=4,fontsize=16)
plt.subplot(322)
plt.plot(np.linspace(-20,240,T),V3q[::100],'b',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-20,240,T),V4q[::100],'r',linewidth=3,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,160)
plt.text(20,20,r'$d=$'+str(lbl[1]),size=25)
plt.legend(loc=4,fontsize=16)
plt.subplot(323)
plt.plot(np.linspace(-20,240,T),V5q[::100],'b',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-20,240,T),V6q[::100],'r',linewidth=3,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,160)
plt.text(20,20,r'$d=$'+str(lbl[2]),size=25)
plt.legend(loc=4,fontsize=16)
plt.subplot(324)
plt.plot(np.linspace(-20,240,T),V7q[::100],'b',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-20,240,T),V8q[::100],'r',linewidth=3,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,160)
plt.text(20,20,r'$d=$'+str(lbl[3]),size=25)
plt.legend(loc=4,fontsize=16)
plt.subplot(325)
plt.plot(np.linspace(-20,240,T),V9q[::100],'b',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-20,240,T),V10q[::100],'r',linewidth=3,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,160)
plt.text(20,20,r'$d=$'+str(lbl[4]),size=25)
plt.legend(loc=4,fontsize=16)
#plt.savefig('EffectOfDistanceC2.eps',dpi=1000,bbox_inches='tight')
<matplotlib.legend.Legend at 0x27c8380d370>
We now test the effect of molar volume for larger species, which permeability is very poor. We extrapolate $v_{Cn}=0.029+0.012(n-1) L.mol^{-1}$.
We observe that the effect of molar volume is very small.
#we vary v_C2 in place of v_Cn
vv1=0.018,0.029,0.041,0.053,0.065 #C2
vv2=0.018,0.029,0.065,0.053,0.065 #C4
vv3=0.018,0.029,0.089,0.053,0.065 #C6
vv4=0.018,0.029,0.161,0.053,0.065 #C12
V0=80
Vc=120
P0=600/20 #illustrated for 14 um
T=260
C10,C20,C1c,C2c=0.0,0.8,0.0,0.08
krel=0 #no reaction
kk=0 #no reaction
#Transport coefficients
P1,P2=0.001*P0,0.001*P0
[V1q,V2q,_,_]=ODEEXCHP5(V0,Vc,vv1,P0,P1,P2,kk,krel,T,C10,C20,C1c,C2c,1)
[V3q,V4q,_,_]=ODEEXCHP5(V0,Vc,vv2,P0,P1,P2,kk,krel,T,C10,C20,C1c,C2c,1)
[V5q,V6q,_,_]=ODEEXCHP5(V0,Vc,vv3,P0,P1,P2,kk,krel,T,C10,C20,C1c,C2c,1)
[V7q,V8q,_,_]=ODEEXCHP5(V0,Vc,vv4,P0,P1,P2,kk,krel,T,C10,C20,C1c,C2c,1)
lbl=[vv1[2],vv2[2],vv3[2],vv4[2]]
matplotlib.pyplot.subplots_adjust(left=None, bottom=None, right=1.9, top=2, wspace=0.3, hspace=0.5)
plt.subplot(221)
plt.plot(np.linspace(-20,240,T),V1q[::100],'b',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-20,240,T),V2q[::100],'r',linewidth=3,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,160)
plt.text(20,20,r'$C_3$',size=25)
plt.legend(loc=4,fontsize=16)
plt.subplot(222)
plt.plot(np.linspace(-20,240,T),V3q[::100],'b',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-20,240,T),V4q[::100],'r',linewidth=3,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,160)
plt.text(20,20,r'$C_4$',size=25)
plt.legend(loc=4,fontsize=16)
plt.subplot(223)
plt.plot(np.linspace(-20,240,T),V5q[::100],'b',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-20,240,T),V6q[::100],'r',linewidth=3,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,160)
plt.text(20,20,r'$C6$',size=25)
plt.legend(loc=4,fontsize=16)
plt.subplot(224)
plt.plot(np.linspace(-20,240,T),V7q[::100],'b',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-20,240,T),V8q[::100],'r',linewidth=3,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,160)
plt.text(20,20,r'$C_{12}$',size=25)
plt.legend(loc=4,fontsize=16)
#plt.savefig('EffectOfMolarVolume.eps',dpi=1000,bbox_inches='tight')
<matplotlib.legend.Legend at 0x27c80a03e50>
First, let us check that when the permeabilities are set to 0, the kinetics of the toy model are recovered.
ratioP2P0 = 0.05
vv=0.018,0.031,0.066,0.099,0.132 #0.018,0.029,0.041,0.053,0.065
V0,Vc=52,75
P0,P1,P2=0,0,0
T=260
Tx=T*10
C10,C20,C1c,C2c=1.0,0.2,1.0,0.0
kk=0.06 #from fit of Nash assay
krel=47*55 #from fit of Nash assay (perfectly agrees with Socha1981), includes Ntot=55 moles/L as a factor because equations expressed with total moles
[VAA,VBB,NNA,NNB]=ODEEXCHP5(V0,Vc,vv,P0,P1,P2,kk,krel,T,C10,C20,C1c,C2c,1)
plt.plot(np.linspace(-20,240,Tx),NNB[1][::10],label='1')
plt.plot(np.linspace(-20,240,Tx),NNB[2][::10],label='2')
plt.plot(np.linspace(-20,240,Tx),NNB[3][::10],label='3')
plt.plot(np.linspace(-20,240,Tx),NNB[4][::10],label='4')
plt.xlim(-20,5)
plt.ylabel('N',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,60)
plt.legend(loc=4,ncol=2,fontsize=16)
<matplotlib.legend.Legend at 0x27c83d2e3a0>
#control: no C2 in droplet B and no exchange, so nothing happens in B
plt.plot(np.linspace(-20,240,Tx),NNA[1][::10],label='1')
plt.plot(np.linspace(-20,240,Tx),NNA[2][::10],label='2')
plt.plot(np.linspace(-20,240,Tx),NNA[3][::10],label='3')
plt.plot(np.linspace(-20,240,Tx),NNA[4][::10],label='4')
plt.xlim(-20,5)
plt.ylabel('N',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.legend(loc=4,ncol=2,fontsize=16)
<matplotlib.legend.Legend at 0x27c851d0d60>
#control volumes constant as no exchange
plt.plot(np.linspace(-20,240,T),VAA[::100],label='1')
plt.plot(np.linspace(-20,240,T),VBB[::100],label='2')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(-0,80)
plt.legend(loc=4,ncol=2,fontsize=16)
<matplotlib.legend.Legend at 0x27c8522b1f0>
We now explore the model with simultaneous exchange and the formose reaction. Note that, unlike in the sugar exchange experiments, the system here contains other species such as methanol (~0.15 M), TMG (1 M), CaCl2 (0.02M).
TMG occupies a considerable volume and is a strong chaotrope. Therefore, the hydrogen bonding in the aqueous droplet phase is destabilized and it becomes energetically less favorable for $H2O$ and hydrophilic species to reside in the aqueous droplet phase.
Since permeability is proportional to the partition constant, our null hypothesis should correspondingly be that H2O, C1 and C2 become more permeable. Below, we retain their original ratios but vary their absolute magnitude to see at what point they produce exchange comparable to what is observed for 20 $\mu m$ separation.
#experimental data formose exchange
import pandas
df = pandas.read_excel(r'C:\Users\Philippe\Desktop\projet_formose\dataDroplets\ED_FR_CTRL.xlsx')
odd30um = df.iloc[0:225,0]
eve30um = df.iloc[0:225,1]
odd20um = df.iloc[0:239,2]
eve20um = df.iloc[0:239,3]
odd14um = df.iloc[0:239,4]
eve14um = df.iloc[0:239,5]
odd8um = df.iloc[0:239,6]
eve8um = df.iloc[0:239,7]
vmean = []
vmean.append(np.mean(odd8um[100:]+eve8um[100:])/2)
vmean.append(np.mean(odd14um[100:]+eve14um[100:])/2)
vmean.append(np.mean(odd20um[100:]+eve20um[100:])/2)
vmean.append(np.mean(odd30um[100:]+eve30um[100:])/2)
display(vmean)
[62.255709852116716, 63.64429393884889, 66.00582562949643, 59.60517896439999]
#compute formose model for all distances,
#and for 2 different values of krel=kadd/kcleave within the uncertainty range observed when fitting Nash bulk formose
distances = [8,14,20,30]
P0list= [4000/x for x in distances]
ratioP2P0 = 0.05
vv=0.018,0.029,0.041,0.053,0.065
V0,Vc=40,80
T=260
Tx=T*10
C10,C20,C1c,C2c=1,0.2,1,0.0
kk=0.06
krel=47*55
Vlistq=[]
for k in range(4):
P0=P0list[k]
P1,P2=P0,ratioP2P0*P0
Vlistq.append(ODEEXCHP5(V0+vmean[k]-60,Vc+vmean[k]-60,vv,P0,P1,P2,kk,krel,T,C10,C20,C1c,C2c,1))
krel2=188*55
Vlistq2=[]
for k in range(4):
P0=P0list[k]
P1,P2=P0,ratioP2P0*P0
Vlistq2.append(ODEEXCHP5(V0+vmean[k]-60,Vc+vmean[k]-60,vv,P0,P1,P2,kk,krel2,T,C10,C20,C1c,C2c,1))
delay=20 #experimental delay of 20 min between the preparation of the droplets and the actual microscopy observation
matplotlib.pyplot.subplots_adjust(left=None, bottom=None, right=1.9, top=2, wspace=0.3, hspace=0.5)
plt.subplot(221)
Time=np.linspace(0,len(odd8um),num=len(odd8um))
plt.plot(Time,odd8um,'ro',alpha=0.3,label=r'data')
plt.plot(Time,eve8um,'bo',alpha=0.3,label=r'data')
k=0
plt.plot(np.linspace(-delay,240,T),Vlistq[k][0][::100],'r',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-delay,240,T),Vlistq[k][1][::100],'b',linewidth=3,label=r'$V^{II}$')
plt.plot(np.linspace(-delay,240,T),Vlistq2[k][0][::100],'r--',linewidth=2,label=r'$V^I$')
plt.plot(np.linspace(-delay,240,T),Vlistq2[k][1][::100],'b--',linewidth=2,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,110)
plt.text(20,20,r'$L=$'+str(distances[k])+' $\mu$m',size=25)
plt.subplot(222)
Time=np.linspace(0,len(odd14um),num=len(odd14um))
plt.plot(Time,odd14um,'ro',alpha=0.3,label=r'data')
plt.plot(Time,eve14um,'bo',alpha=0.3,label=r'data')
k=1
plt.plot(np.linspace(-delay,240,T),Vlistq[k][0][::100],'r',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-delay,240,T),Vlistq[k][1][::100],'b',linewidth=3,label=r'$V^{II}$')
plt.plot(np.linspace(-delay,240,T),Vlistq2[k][0][::100],'r--',linewidth=2,label=r'$V^I$')
plt.plot(np.linspace(-delay,240,T),Vlistq2[k][1][::100],'b--',linewidth=2,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,110)
plt.text(20,20,r'$L=$'+str(distances[k])+' $\mu$m',size=25)
plt.subplot(223)
Time=np.linspace(0,len(odd20um),num=len(odd20um))
plt.plot(Time,odd20um,'ro',alpha=0.3,label=r'data')
plt.plot(Time,eve20um,'bo',alpha=0.3,label=r'data')
k=2
plt.plot(np.linspace(-delay,240,T),Vlistq[k][0][::100],'r',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-delay,240,T),Vlistq[k][1][::100],'b',linewidth=3,label=r'$V^{II}$')
plt.plot(np.linspace(-delay,240,T),Vlistq2[k][0][::100],'r--',linewidth=2,label=r'$V^I$')
plt.plot(np.linspace(-delay,240,T),Vlistq2[k][1][::100],'b--',linewidth=2,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,110)
plt.text(20,20,r'$L=$'+str(distances[k])+' $\mu$m',size=25)
plt.subplot(224)
Time=np.linspace(0,len(odd30um),num=len(odd30um))
plt.plot(Time,odd30um,'ro',alpha=0.3,label=r'data')
plt.plot(Time,eve30um,'bo',alpha=0.3,label=r'data')
k=3
plt.plot(np.linspace(-delay,240,T),Vlistq[k][0][::100],'r',linewidth=3,label=r'$V^I$')
plt.plot(np.linspace(-delay,240,T),Vlistq[k][1][::100],'b',linewidth=3,label=r'$V^{II}$')
plt.plot(np.linspace(-delay,240,T),Vlistq2[k][0][::100],'r--',linewidth=2,label=r'$V^I$')
plt.plot(np.linspace(-delay,240,T),Vlistq2[k][1][::100],'b--',linewidth=2,label=r'$V^{II}$')
plt.ylabel('Volume',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
plt.ylim(0,110)
plt.text(20,20,r'$L=$'+str(distances[k])+' $\mu$m',size=25)
#plt.savefig('FormoseDroplets.eps',dpi=1000,bbox_inches='tight')
Text(20, 20, '$L=$30 $\\mu$m')
We now examine the dynamics of species in each of the coupled droplets.
V0,Vc=42,80
ratioP2P0 = 0.05
vv=0.018,0.029,0.041,0.053,0.065
P0=50.0 #12.5
T=260
Tx=T*10
C10,C20,C1c,C2c=1.0,0.2,1.0,0.0
kk=0.06
krel=55.0*47 #PN changed from 55.0*187.5
P1,P2=1.0*P0,ratioP2P0*P0
[VAA,VBB,NNA,NNB]=ODEEXCHP5(V0,Vc,vv,P0,P1,P2,kk,krel,T,C10,C20,C1c,C2c,1)
matplotlib.pyplot.subplots_adjust(left=None, bottom=None, right=2.5, top=2, wspace=0.3, hspace=0.5)
plt.subplot(211)
plt.plot(np.linspace(-20,240,Tx),NNA[1][::10],label='1')
plt.plot(np.linspace(-20,240,Tx),NNA[2][::10],label='2')
plt.plot(np.linspace(-20,240,Tx),NNA[3][::10],label='3')
plt.plot(np.linspace(-20,240,Tx),NNA[4][::10],label='4')
#plt.xlim(-20,5)
plt.ylabel('N',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
#plt.ylim(0,int(1.1*max([max(Vv3),max(Vv4)])))
plt.legend(loc=4,ncol=2,fontsize=16)
plt.subplot(212)
plt.plot(np.linspace(-20,240,Tx),NNB[1][::10],label='1')
plt.plot(np.linspace(-20,240,Tx),NNB[2][::10],label='2')
plt.plot(np.linspace(-20,240,Tx),NNB[3][::10],label='3')
plt.plot(np.linspace(-20,240,Tx),NNB[4][::10],label='4')
#plt.xlim(-20,5)
plt.ylabel('N',size=20)
plt.xlabel('Time (min)',size=20)
plt.xticks(size=16)
plt.yticks(size=16)
#plt.ylim(0,int(1.1*max([max(Vv3),max(Vv4)])))
plt.legend(loc=4,ncol=2,fontsize=16)
<matplotlib.legend.Legend at 0x27c80a976a0>
Let us now explore more systematically various regimes afforded by varying the timescales of the reaction and transport.
#difference with ODEEXCHP5: outputs only the final volume of droplet A, which saves time to compute the phase diagram
def ODEEXVF(VA,VB,vv,PH2O,P1,P2,kk,krel,T,C1A,C2A,C1B,C2B,nplus):
# ODE solver parameters
TT=1
abserr = 1.0e-10
relerr = 1.0e-11
stoptime = T*1.0
numpoints = T*100*nplus
# Create the time samples for the output of the ODE solver.
t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]
k12=kk*krel #C1+C2 aldol rate
k13=kk*krel #C1+C3 aldol rate
k4=kk #C4 cleavage by retro-aldol
k22=kk*krel #C2+C2->C4 aldol
Ninert=0.00
[vvH2O,vvC1,vvC2,vvC3,vvC4]=vv
k=[P1,P2,k12,k13,k4,PH2O,k22,Ninert,vvH2O,vvC1,vvC2,vvC3,vvC4]
NL0=zeros(10)
NL0[0]=VA*(1.0-vvC1*C1A-vvC2*C2A)/vvH2O #NH2O in A
NL0[1]=VA*C1A #NC1 in A
NL0[2]=VA*C2A #NC2 in A
NL0[3]=0.0 #NC3 in A
NL0[4]=0.0 #NC4 in A
NL0[5]=VB *(1.0-C1B*vvC1-C2B*vvC2)/vvH2O #NH2O in B
NL0[6]=VB * C1B #NC1 in B
NL0[7]=VB * C2B #NC1 in B
NL0[8]=0.0 #NC3 in B
NL0[9]=0.0 #NC4 in B
# Call the ODE solver.
NLF = odeint(vexch7, NL0, t, args=(k,),atol=abserr, rtol=relerr)
NFin=numpoints-1
VAF=(NLF[NFin][0]*vvH2O + NLF[NFin][1]*vvC1 + NLF[NFin][2]*vvC2 + NLF[NFin][3]*vvC3 + NLF[NFin][4]*vvC4)
return VAF
NM=51
PXL=np.logspace(-2,4,NM)#ratios P1/P2, will be x axis in phase diagram
kkL=np.logspace(-4,3,NM)#typical reaction rate, will be y axis in phase diagram
NPL=np.concatenate((np.linspace(1,1,NM-10),np.linspace(10,10,5),np.linspace(50,50,5)),axis=0)
P0=200
P2=0.05*P0
vv=0.018,0.029,0.041,0.053,0.065
VA0,VB0=40,80
T=260
Tx=T*10
C1A,C2A,C1B,C2B=1.0,0.2,1.0,0.0
kk=0.06
krel=55*47
#Computes volume fold change of droplet A containing some C2 at the beginning
DVLt=[]
#stored along 'columns' for varying PX=P1/P2=TC2/TC1 (the array within the array)
#stored along 'rows' dimension for varying kk
for i in range(NM):
kk=kkL[i] #rate of the limiting step of the reaction C4->2C2, in min^-1
DVLt.append([]) #kk increases along rows of DVLt
for PX in PXL:
P1=PX*P2
nplus=int(NPL[i]) #extra points, needed for integration when rates are high
VAF=ODEEXVF(VA0,VB0,vv,P0,P1,P2,kk,krel,T,C1A,C2A,C1B,C2B,nplus)
DVLt[i].append(VAF)#Final over average initial volume
To obtain a characteristic exchange time of C2, we consider a droplet with C2 at concentration $c$ coupled to a droplet without C2. The flux is $P_{C2} A_{red} x_{C2}$ in $mol.min^{-1}$ with $x_{C2}=n_{C2}/n_{tot} \approx n_{C2}/n_{H2O}$ and transports a characteristic amount of C2 of $n_{C2}$ in moles. In fine, $\tau_{C2} \approx c_{H2O} V / (P_{C2} A_{red})\approx c_{H2O}D/P_{C2}$ where $D$ is the droplet diameter.
area = VA0**(2.0/3.0)*VB0**(2.0/3.0)/(VA0**(2.0/3.0)+VB0**(2.0/3.0))
tauC2=55*VA0**(1/3)/P2 #characteristic time of C2 exchange in minutes
display(log10(tauC2))
display(tauC2)
1.2743826866035646
18.809735413443665
y = [log10(u*tauC2) for u in kkL]
X, Y = np.meshgrid(log10(PXL),y)
Z=[[x/VA0 for x in lst] for lst in DVLt]
levels = np.arange(0.8, 2.4, 0.2)
fig, ax = plt.subplots()
fig.set_size_inches(5,4)
CS = ax.contour(X, Y, Z,levels)
plt.xlabel(r'$P_{C_1} / P_{C_2}$',size=16)
plt.ylabel(r'$k_{cleave} c_0 D / P_{C_2}$',size=16)
ax.clabel(CS, inline=True, fontsize=9,fmt='%1.1f')
xexp = 1/0.05
yexp = 0.06*tauC2
plt.plot(log10(xexp),log10(yexp),'ko')
#plt.savefig('PhaseDiagram.eps',dpi=1000,bbox_inches='tight')
[<matplotlib.lines.Line2D at 0x27c85c82760>]
#kcleave*tauC2
display(0.06*tauC2)
1.12858412480662