Multi-tube analysis example

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

In [1]:
# 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
Out[1]:
Complex results:
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
Concentration results:
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
In [2]:
# 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)
In [3]:
# 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)
In [4]:
# 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)
In [5]:
# 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)