Analyze multiple test tube ensembles that model different states of HCR polymerization (Dirks and Pierce, Proc Natl Acad Sci USA, 2004)
Strand species:
Initiator i1
Hairpin h1
Hairpin h2
Tube t1:
i1: 1e-9 M
h1: 1e-8 M
h2: 1e-8 M
All complexes of up to 3 strands plus the cognate 4-mer (polymer ending with h1)
Tube t2:
i1: 1e-9 M
h1: 1e-8 M
h2: 1e-8 M
All complexes of up to 3 strands plus the cognate 4-mer and 5-mer (polymer ending with h2)
Material: DNA
Temperature: 23 C
# Import NUPACK Python module
from nupack import *
# Define physical model
my_model = Model(material='dna', celsius=23)
# Define strand species
i1 = Strand('AGTCTAGGATTCGGCGTGGGTTAA', name='i1')
h1 = Strand('TTAACCCACGCCGAATCCTAGACTCAAAGTAGTCTAGGATTCGGCGTG', name='h1')
h2 = Strand('AGTCTAGGATTCGGCGTGGGTTAACACGCCGAATCCTAGACTACTTTG', name='h2')
# Define two tube ensembles containing strands at specified concentrations interacting
# to form all complexes of up to 3 strands plus additional specified complexes
t1 = Tube(strands={i1: 1e-9, h1: 1e-8, h2: 1e-8},
complexes=SetSpec(max_size=3,
include=[[i1, h1, h1, h2]]),
name='tube t1')
t2 = Tube(strands={i1: 1e-9, h1: 1e-8, h2: 1e-8},
complexes=SetSpec(max_size=3,
include=[[i1, h1, h1, h2],
[i1, h1, h1, h2, h2]]),
name='tube t2')
# Analyze the tube ensembles
# Calculate pfunc (default), pairs, mfe for each complex
# Calculate equilibrium complex concentrations (default) for each tube ensemble
# Since pairs is specified, calculate ensemble pair fractions for each tube ensemble
tube_result = tube_analysis(tubes=[t1, t2], compute=['pairs', 'mfe'], model=my_model)
tube_result
Complex | Pfunc | ΔG (kcal/mol) | MFE (kcal/mol) |
---|---|---|---|
(h1) | 3.8928e+20 | -27.901 | -27.740 |
(h2) | 7.5022e+20 | -28.287 | -27.979 |
(i1) | 6.0965e+0 | -1.064 | 0.000 |
(h1+h1) | 2.8860e+44 | -60.247 | -59.602 |
(h2+h1) | 8.9020e+44 | -60.910 | -59.538 |
(h2+h2) | 1.1453e+45 | -61.058 | -59.478 |
(i1+h1) | 6.1982e+31 | -43.081 | -41.379 |
(i1+h2) | 1.9312e+26 | -35.619 | -33.284 |
(i1+i1) | 1.1115e+9 | -12.258 | -11.507 |
(h1+h1+h1) | 1.1499e+68 | -92.227 | -90.688 |
(h2+h1+h1) | 9.8849e+68 | -93.493 | -91.242 |
(h2+h2+h1) | 3.6082e+69 | -94.255 | -91.951 |
(h2+h2+h2) | 2.4956e+69 | -94.038 | -91.454 |
(i1+h1+h1) | 3.9049e+57 | -78.041 | -76.006 |
(i1+h1+h2) | 7.8258e+62 | -85.225 | -83.629 |
(i1+h2+h1) | 6.8345e+54 | -74.305 | -71.222 |
(i1+h2+h2) | 3.2841e+51 | -69.808 | -66.297 |
(i1+i1+h1) | 4.0209e+40 | -55.022 | -52.886 |
(i1+i1+h2) | 6.0106e+34 | -47.128 | -44.791 |
(i1+i1+i1) | 2.5831e+15 | -20.885 | -19.339 |
(i1+h1+h1+h2) | 2.2947e+94 | -127.866 | -126.035 |
(i1+h1+h1+h2+h2) | 2.8973e+125 | -170.010 | -168.285 |
Complex | tube t1 (M) | Complex | tube t2 (M) | ||
---|---|---|---|---|---|
(h2) | 9.037e-09 | (h2) | 8.323e-09 | ||
(h1) | 8.124e-09 | (h1) | 8.040e-09 | ||
(i1+h1+h1+h2) | 8.834e-10 | (i1+h1+h1+h2+h2) | 6.899e-10 | ||
(i1+h1+h2) | 7.993e-11 | (i1+h1+h1+h2) | 2.727e-10 | ||
(i1+h1) | 2.910e-11 | (i1+h1+h2) | 2.494e-11 | ||
(i1) | 7.595e-12 | (i1+h1) | 9.858e-12 | ||
(h2+h1) | 4.041e-15 | (i1) | 2.600e-12 | ||
(h2+h2) | 3.001e-15 | (h2+h1) | 3.683e-15 | ||
(h1+h1) | 2.270e-15 | (h2+h2) | 2.545e-15 | ||
(i1+h1+h1) | 6.910e-16 | (h1+h1) | 2.223e-15 | ||
(i1+i1+h1) | 4.247e-16 | (i1+h1+h1) | 2.317e-16 | ||
(i1+h2) | 5.233e-17 | (i1+i1+h1) | 4.925e-17 | ||
(i1+i1) | 3.115e-17 | (i1+h2) | 1.650e-17 | ||
(i1+h2+h1) | 6.981e-19 | (i1+i1) | 3.650e-18 | ||
(h2+h2+h1) | 3.563e-21 | (i1+h2+h1) | 2.178e-19 | ||
(h2+h1+h1) | 1.691e-21 | (h2+h2+h1) | 2.991e-21 | ||
(h2+h2+h2) | 1.422e-21 | (h2+h1+h1) | 1.526e-21 | ||
(i1+i1+h2) | 3.664e-22 | (h2+h2+h2) | 1.111e-21 | ||
(h1+h1+h1) | 3.409e-22 | (h1+h1+h1) | 3.304e-22 | ||
(i1+h2+h2) | 1.936e-22 | (i1+h2+h2) | 5.621e-23 | ||
(i1+i1+i1) | 1.629e-24 | (i1+i1+h2) | 3.954e-23 | ||
(i1+i1+i1) | 6.531e-26 |
# MFE proxy structure for odd polymer ending in h1
odd_polymer_result = tube_result['(i1+h1+h1+h2)']
print('\nMFE proxy structure(s) for odd polymer ending in h1:')
for i, s in enumerate(odd_polymer_result.mfe):
print(' %2d: %s (%.2f kcal/mol)' % (i, s.structure, s.energy))
# Plot the equilibrium pair probability matrix for odd polymer ending in h1
import matplotlib.pyplot as plt
plt.imshow(odd_polymer_result.pairs.to_array())
plt.xlabel('Base index')
plt.ylabel('Base index')
plt.title('Pair probabilities for complex (i1+h1+h1+h2)')
plt.colorbar()
plt.clim(0, 1)
MFE proxy structure(s) for odd polymer ending in h1: 0: ((((((((((((((((((((((((+))))))))))))))))))))))))((((((((((((((((((((((((+(((((((((((((((((((((((((..((....)).)...........+)))))))))))))))))))))))))))))))))))))))))))))))) (-126.03 kcal/mol)
# MFE proxy structure for even polymer ending in h2
even_polymer_result = tube_result['(i1+h1+h1+h2+h2)']
print('\nMFE proxy structure(s) for even polymer ending in h2:')
for i, s in enumerate(even_polymer_result.mfe):
print(' %2d: %s (%.2f kcal/mol)' % (i, s.structure, s.energy))
# Plot the equilibrium pair probability matrix for even polymer ending in h2
import matplotlib.pyplot as plt
plt.imshow(even_polymer_result.pairs.to_array())
plt.xlabel('Base index')
plt.ylabel('Base index')
plt.title('Pair probabilities for complex (i1+h1+h1+h2+h2)')
plt.colorbar()
plt.clim(0, 1)
MFE proxy structure(s) for even polymer ending in h2: 0: ((((((((((((((((((((((((+))))))))))))))))))))))))((((((((((((((((((((((((+((((((((((((((((((((((((((((((((((((((((((((((((+........................))))))))))))))))))))))))+)))))))))))))))))))))))))))))))))))))))))))))))) (-168.29 kcal/mol)
# Plot the ensemble pair fractions for tube t1
plt.imshow(tube_result[t1].ensemble_pair_fractions.to_array())
plt.xlabel('Base index')
plt.ylabel('Base index')
plt.title('Ensemble pair fractions for ' + ', '.join(s.name for s in t1.strands))
plt.colorbar()
plt.clim(0, 1)
# Plot the ensemble pair fractions for tube t2
plt.imshow(tube_result[t2].ensemble_pair_fractions.to_array())
plt.xlabel('Base index')
plt.ylabel('Base index')
plt.title('Ensemble pair fractions for ' + ', '.join(s.name for s in t2.strands))
plt.colorbar()
plt.clim(0, 1)