<jupyter_start><jupyter_text>Calculating $R^2$ of kinked WLC for heterogenous linker lengths<jupyter_code>%run init.ipy
%matplotlib inline
#Generate code to randomly sample between 35 and 36 bp linker lengths, 20 unwraps
plt.figure()
for i in range(5):
    links = np.random.randint(35, 37, 10000)
    r2, ldna, m = wlc.kinked_WLC_no_translation(links, 'blah', plotfig=False)
    plt.loglog(ldna, r2)
    plt.xlabel(r'$R_{max}$ [nm]')
    plt.ylabel(r'$\langle{R^2}\rangle$ [nm$^2$]')
    plt.title('Random linkers 35-36 bp')

plt.text(0.10, 0.80, f'Kuhn length: {np.mean(kuhns):.2f} nm', transform=ax.transAxes, size=12)
plt.subplots_adjust(bottom = 0.15, top = 0.9, left = 0.13, right = 0.95)
plt.savefig('RandomLinkers35to36bp.pdf')
r2 = np.zeros((5, 10001))
kuhns = np.zeros(5)
ldna = np.zeros_like(r2)
for i in range(5):
    links = np.random.randint(35, 50, 10000)
    r2d, ldnad, kuhnsd = wlc.kinked_WLC_no_translation(links, 'blah', plotfig=False)
    r2[i, :] = r2d
    ldna[i, :] = ldnad
    kuhns[i] = kuhnsd

print(f'Mean Kuhn length in nm: {np.mean(kuhns):.2f}')<jupyter_output>/Users/deepti/Documents/Spakowitz/nuc_chain/nuc_chain/fluctuations.py:39: RuntimeWarning: invalid value encountered in sqrt
  al  = lambda l, j: np.sqrt((l-j)*(l+j)/(4*l**2 - 1)) #ladder coefficients alpha<jupyter_text>Above, we have 5 different chain realizations with random linkers between 35-49 bp (excluding unwrapping, which by default is set to 10 bp on either side, so total linker length is 55-69 bp). The kinks seem to eliminate the typical 2-slope bend you see in $R^2$ vs. $R_{max}$ curves (see below). The length scales where random linkers matter extends up to 1000 - 10,000 bp, which is a significant number of nucleosomes.<jupyter_code>plt.figure()
ax = plt.subplot(111)
for i in range(5):
    plt.loglog(ldna[i,:], r2[i,:])

plt.xlabel('$R_{max}$ [nm]')
plt.ylabel('$<R^2>$ [nm$^2$]')
plt.text(0.10, 0.80, f'Mean Kuhn length: {np.mean(kuhns):.2f} nm', transform=ax.transAxes, size=12)
plt.title('Random linkers 35-49 bp')
plt.subplots_adjust(bottom = 0.15, top = 0.9, left = 0.13, right = 0.95)
plt.savefig('RandomLinkers35to50bp')
print(f'Mean Kuhn length in nm: {np.mean(kuhns):.2f}')<jupyter_output>Mean Kuhn length in nm: 2.49<jupyter_text>1. Test no kinks / no twist limitsFor a normal WLC with no twist and no kinks, the analytical formula for $$ is given by\begin{equation*}\langle{R^2}\rangle = 2l_pR_{max}(1 - \frac{l_p}{R_{max}}(1 - e^{-R_{max}/l_p})\end{equation*}Notice this formula has no dependence on the twist persistance length. To test against this formula, we set Euler angles for all kinks to be (0, 0, 0), resulting in a normal WLC.Our plots are on a log-log scale. What does this mean?In the ideal chain limit, $R_{max} >> l_p$, so $\langle{R^2}\rangle = 2l_pR_{max}$. Taking the log of both sides, we get $log(\langle{R^2}\rangle) = log(R_{max}) + log(2l_p)$. Thus, in the long chain limit, the slope of the log-log plot will always be 1, and the persistance length can be extracted from the intercept.In the rigid rod limit, $R_{max} << l_p$, so $\langle{R^2}\rangle = R_{max}^2$. But this is for a normal worm-like chain (our kinked chain won't necessarily follow this behavior).<jupyter_code>#Hard code euler angles=0 to test agreement with normal WLC (no kinks, totally straight polymer)
%matplotlib inline

def r2wlc(ldna, lp):
    """Analytical formula for R^2 of WLC."""
    return 2*(lp*ldna - lp**2 + lp**2*np.exp(-ldna/lp))

links = 35*np.ones(100000)
r2nokinks, ldnanokinks, mnokinks = wlc.kinked_WLC_no_translation(links, lt=0.01*wlc.default_lt)
print(mnokinks)
r2theory = r2wlc(ldnanokinks, wlc.default_lp*ncg.dna_params['lpb'])

plt.figure()
ax = plt.subplot(111)
plt.loglog(ldnanokinks, r2theory)
plt.loglog(ldnanokinks, r2nokinks)

plt.legend([r'$\alpha, \gamma, \beta = 0$', 'Analytical WLC'])
plt.xlabel('$R_{max}$ [nm]')
plt.ylabel('$<R^2>$ [nm$^2$]')
plt.text(0.5, 0.2, f'Kuhn length: {mnokinks:.2f} nm', transform=ax.transAxes, size=12)
plt.title('WLC limit: no kinks, no twist')
plt.subplots_adjust(bottom = 0.15, top = 0.9, left = 0.13, right = 0.95)
plt.savefig('R2_no_kinks_35bplinkers.pdf')
#try adding more data points by discretizing by 1 bp, instead of by 35
links = np.ones(100000)
r2tw, ldna, m = wlc.kinked_WLC_no_translation(links, unwraps=0)
r2theory = r2wlc(ldna, wlc.default_lp*ncg.dna_params['lpb'])

print(ldna[0:5])
plt.figure()
ax = plt.subplot(111)
plt.loglog(ldna, r2theory)
plt.loglog(ldna, r2tw)
plt.legend([r'$\alpha, \gamma, \beta = 0$', 'Analytical WLC'])
plt.xlabel('$R_{max}$ [nm]')
plt.ylabel('$<R^2>$ [nm$^2$]')
plt.text(0.5, 0.2, f'Kuhn length: {mnokinks:.2f} nm', transform=ax.transAxes, size=12)
plt.title('WLC limit: no kinks, no twist')
plt.subplots_adjust(bottom = 0.15, top = 0.9, left = 0.13, right = 0.95)
plt.savefig('R2_no_kinks_1bplinkers.pdf')<jupyter_output>[0.332 0.664 0.996 1.328 1.66 ]<jupyter_text>When we put in 1 bp long linkers with Euler angles set to 0 for the kink, we find that as expected, the M matrices are just the identity. I also printed out the g matrices and they match my analytical expectations.<jupyter_code>links = np.ones(100000)
#crank down the twist persistence length 
r2p, ldna, m = wlc.kinked_WLC_no_translation(links, 'R2_no_kinks_Lt0.001_1bp_linkers.pdf', lt=0.001)
print(m)
r2theory = r2wlc(ldna, wlc.default_lp)
plt.loglog(ldna, r2theory)
plt.legend(['WLC no kinks, lt=0.001', 'Analytical WLC'])
plt.xlabel('Length of chain in basepairs')
plt.ylabel(r'$<R^2>$')<jupyter_output>/Users/deepti/Documents/Spakowitz/nuc_chain/nuc_chain/fluctuations.py:39: RuntimeWarning: invalid value encountered in sqrt
  al  = lambda l, j: np.sqrt((l-j)*(l+j)/(4*l**2 - 1)) #ladder coefficients alpha<jupyter_text>Cranking down the twist persistance length doesn't change anything. This is to be expected because there is no twist-bend coupling. Without kinks, twist doesn't change $R^2$, which only depends on the tangent vector and not the normal vectors (which are the only ones that get rotated by twist, and therefore the only ones that are affected by the twist persistance length).<jupyter_code>links = np.ones(100000)
#crank down the twist persistence length 
r2p, ldna, m = wlc.kinked_WLC_no_translation(links, 'R2_no_kinks_10Lt_1bp_linkers.pdf', lt=10*wlc.default_lt)
print(m)
r2theory = r2wlc(ldna, wlc.default_lp)
plt.loglog(ldna, r2theory)
plt.legend(['WLC no kinks, 10Lt', 'Analytical WLC'])
plt.xlabel('Length of chain in basepairs')
plt.ylabel(r'$<R^2>$')<jupyter_output>/Users/deepti/Documents/Spakowitz/nuc_chain/nuc_chain/fluctuations.py:39: RuntimeWarning: invalid value encountered in sqrt
  al  = lambda l, j: np.sqrt((l-j)*(l+j)/(4*l**2 - 1)) #ladder coefficients alpha<jupyter_text>Neither does cranking up the twist persistence length...  2. Test rigid rod limit.In the rigid rod limit, $l_p$ and $l_t$ go to infinity, so only term in Hamiltonian is $-i\tau\frac{\partial}{\partial \psi}$. In this limit (no fluctuations), $\langle{R^2}\rangle$ should match the pure geometric picture.<jupyter_code>links = 35*np.ones(100000)
r2p, ldnap, mp = wlc.kinked_WLC_no_translation(links, 'R2_kinks_10Lt_35bp_linkers.pdf', plotfig=False, lt=10*wlc.default_lt)
r2t, ldnat, mt = wlc.kinked_WLC_no_translation(links, 'R2_kinks_0.01Lt_35bp_linkers.pdf', plotfig=False, lt=0.01*wlc.default_lt)
print(f'Kuhn length in bp for 10Lt: {mp:.2f}')
print(f'Kuhn length in bp for .01Lt: {mt:.2f}')

plt.figure()
ax = plt.subplot(111)
plt.loglog(ldnap, r2p)
plt.loglog(ldnat, r2t)
plt.legend(['10$l_t$', '0.01$l_t$'])
plt.xlabel('$R_{max}$ [nm]')
plt.ylabel('$<R^2>$ [nm$^2$]')
plt.title('Effect of twist persistance length for 35bp linkers')
plt.subplots_adjust(bottom = 0.15, top = 0.9, left = 0.13, right = 0.95)
plt.savefig('R2_kinks_35bplinkers_10lt_0.01lt.pdf')<jupyter_output>/Users/deepti/Documents/Spakowitz/nuc_chain/nuc_chain/fluctuations.py:39: RuntimeWarning: invalid value encountered in sqrt
  al  = lambda l, j: np.sqrt((l-j)*(l+j)/(4*l**2 - 1)) #ladder coefficients alpha<jupyter_text>The above plot confirms that changing the twist persistance length DOES matter for the kinked chain. Here, default $l_t$ is set to be 100 nm, whereas default $l_p$ is 50 nm.<jupyter_code>#Now let's generate R^2 for 
links = 35*np.ones(100000)
rots, pos = ncg.minimum_energy_no_sterics_linker_only(links)
#calculate R^2 in nm
r2 = np.array([np.linalg.norm(pos[i, :] - pos[0, :])**2 for i in range(pos.shape[0])])
ldna = np.array([i*55 for i in range(len(r2))]) * ncg.dna_params['lpb'] 

plt.loglog(ldna, r2)
plt.xlabel(r'$R_{max}$ [nm]')
plt.ylabel(r'$\langle{R^2}\rangle$ [nm$^2$]')
#compare against fluctuation code with high persistance length and high twist persistance length
plt.figure()
plt.loglog(ldna, np.array(r2))
plt.xlabel('$R_{max}$ [nm]')
plt.ylabel('$<R^2>$ [nm$^2$]')

r2rr, ldnarr, m = wlc.kinked_WLC_no_translation(links, 'blah', plotfig=False, lp=1000*wlc.default_lp,
                                             lt=1000*wlc.default_lt)
print(f'Kuhn length in nm: {m:.2f}')
plt.loglog(ldnarr, r2rr)
plt.legend(['geometric picture', '$1000l_p$, $1000l_t$'])
plt.title('Kinked WLC in rigid rod limit')
plt.subplots_adjust(bottom = 0.15, top = 0.9, left = 0.13, right = 0.95)
plt.savefig('R2_rigid_rod_limit_35bplinkers.pdf')
ldnarr[1000]
#seems to be off. try re-scaling twist density by 0.9
links = 35*np.ones(100000)
rots, pos = ncg.minimum_energy_no_sterics_linker_only(links, tau_d=ncg.dna_params['tau_d'])
r2 = [np.linalg.norm(pos[i, :] - pos[0, :])**2 for i in range(pos.shape[0])]
ldna = [i*55 for i in range(len(r2))]
plt.loglog(ldna, r2)
plt.xlabel(r'$R_{max}$ in bp')
plt.ylabel(r'$\langle{R^2}\rangle$ in bp$^2$')
plt.loglog(ldnarr, r2rr)
plt.legend(['geometric picture', '$1000l_p$, $1000l_t$'])
#in linear space
plt.plot(ldna, r2)
plt.plot(ldnarr, r2rr)
plt.figure(figsize=(10,6))
plt.loglog(ldna, 9*np.array(r2))
plt.loglog(ldnarr, 1*r2rr)<jupyter_output><empty_output><jupyter_text>Effect of Heterogenous Linkers on Kuhn lengthNext, we compare our heterogenous kinked chain to kinked chains with fixed linkers.<jupyter_code>#First calculate Kuhn length as a function of linker length
links_35to50 = np.arange(35, 50)
kuhns_35to50 = np.zeros_like(links_35to50)
for i in range(len(links_35to50)):
    linkers = np.tile(links_35to50[i], 100000)
    r2, ldna, kuhn = wlc.kinked_WLC_no_translation(linkers, 'blah', plotfig=False)
    kuhns_35to50[i] = kuhn

mean_kuhn_35to50 = np.tile(np.mean(kuhns_35to50), len(links_35to50))
plt.plot(links_35to50, mean_kuhn_35to50, 'r--')
plt.plot(links_35to50, kuhns_35to50, '-o')
plt.xlabel('Linker length [bp]')
plt.ylabel('Kuhn length [nm]')
plt.title('Kuhn length for kinked WLCs with fixed linkers')
plt.savefig('kuhn_length_vs_linkers35to49bp_default_unwrapping.pdf')
#First calculate Kuhn length as a function of linker length with no unwrapping
links = np.arange(1, 251)
kuhns = np.zeros_like(links)
for i in range(len(links)):
    linkers = np.tile(links[i], 100000)
    r2, ldna, kuhn = wlc.kinked_WLC_no_translation(linkers, 'blah', plotfig=False, unwraps=0)
    kuhns[i] = kuhn

plt.plot(links, kuhns, '-o', markersize=4, linewidth=1)
plt.xlabel('Linker length [bp]')
plt.ylabel('Kuhn length [nm]')
plt.title('Kuhn length for kinked WLCs with fixed linkers')
plt.savefig('kuhn_length_vs_linkers1to250bp_no_unwrapping.pdf')
#First calculate Kuhn length as a function of linker length with default unwrapping (10 bp on each side)
links_wrap = np.arange(1, 251)
kuhns_wrap = np.zeros_like(links_wrap)
for i in range(len(links_wrap)):
    linkers = np.tile(links_wrap[i], 100000)
    r2, ldna, kuhn = wlc.kinked_WLC_no_translation(linkers, 'blah', plotfig=False)
    kuhns_wrap[i] = kuhn

plt.plot(links_wrap, kuhns_wrap, '-o', markersize=4, linewidth=1)
plt.xlabel('Linker length [bp]')
plt.ylabel('Kuhn length [nm]')
plt.title('Kuhn length for kinked WLCs with fixed linkers')
plt.savefig('kuhn_length_vs_linkers1to250bp_default_unwrapping.pdf')<jupyter_output>/Users/deepti/Documents/Spakowitz/nuc_chain/nuc_chain/fluctuations.py:39: RuntimeWarning: invalid value encountered in sqrt
  al  = lambda l, j: np.sqrt((l-j)*(l+j)/(4*l**2 - 1)) #ladder coefficients alpha<jupyter_text>As expected, the Kuhn length is 10.5 bp periodic with linker length. The average Kuhn length increases linearly with linker length (just as it does in the freely jointed case). In addition as the linkers get longer, the variation in Kuhn length due to the 10.5 bp-periodic kinks reduces (which makes sense).<jupyter_code>#first check if index 5000 is a fair number to use for this wide range of linker lengths
linkers = np.tile(250, 100000)
r2, ldna, kuhn = wlc.kinked_WLC_no_translation(linkers, 'R2_250bp_linkers.pdf', plotfig=True, unwraps=145)
print(ldna[0:5]/ncg.dna_params['lpb'])
print(kuhn)
plt.plot(ldna[5000:], r2[5000:])<jupyter_output>[ 395.  790. 1185. 1580. 1975.]
99.89389432611321<jupyter_text>Decoding relationship between fixed linker chains and heterogenous chains__Hypothesis 1:__According to "Exact theory of kinkable elastic polymers" by Wiggins and Nelson,\begin{equation*}\frac{1}{\xi_{eff}} = \frac{1}{\xi} + \frac{1}{\xi_{kink}}\end{equation*}where $\xi_{kink} = \zeta^{-1}$, and $\zeta$ is the kink density. I tried this formula using $\zeta$ =1/(linker length) but it didn't match the Kuhn lengths I was getting for chains with fixed linkers. This is because Wiggins/Nelson's model included an energy associated with the kink, whereas in our case the kink is determined.__Hypothesis 2: Andy's guess__\begin{equation*}\frac{1}{\xi_{eff}} = \langle{\frac{1}{\xi}}\rangle = \frac{1}{N}\sum{\frac{1}{\xi_i}}\end{equation*}where N is the number of unique linker lengths and $\xi_i$ is the persistance length for the kinked chain corresponding to each of the fixed linker lengths. Turns out this formula works pretty darn well!<jupyter_code>#determine kink densities for linkers 35 to 50 bp with default wrapping
#kink density is number of kinks per unit length (aka 1/linker length)

links_35to50_nm = links_35to50 * ncg.dna_params['lpb']
eff_kuhns_35to50 = 2*((1/50.) + (1/links_35to50_nm))**(-1)
print(eff_kuhns_35to50)
print(kuhns_35to50)
#ahh this sort of matches Andy's formula!!!
print(np.mean(1/kuhns_35to50)**(-1))
#doesn't work for 35 vs 36 bp though
print(np.mean(1/kuhns_35to50[0:2])**(-1))
mean_kuhn_35to50_no_unwrap = np.tile(np.mean(kuhns[35:47]), len(links[35:47]))
plt.plot(links[35:47], mean_kuhn_35to50_no_unwrap, 'r--')
plt.plot(links[35:47], kuhns[35:47], '-o')
plt.xlabel('Linker length [bp]')
plt.ylabel('Kuhn length [nm]')
plt.title('Kuhn length for kinked WLCs 36-47 bp linkers')
plt.savefig('kuhn_length_vs_linkers36to47bp_no_unwrapping.pdf')
r2_36to47_no_unwrap = np.zeros((10, 10000))
kuhns_36to47_no_unwrap = np.zeros(10)
ldna_36to47_no_unwrap = np.zeros_like(r2_36to47_no_unwrap)
for i in range(10):
    links = np.random.randint(36, 48, 10000)
    r2d, ldnad, kuhnsd = wlc.kinked_WLC_no_translation(links, 'blah', plotfig=False, unwraps=0)
    r2_36to47_no_unwrap[i, :] = r2d
    ldna_36to47_no_unwrap[i, :] = ldnad
    kuhns_36to47_no_unwrap[i] = kuhnsd

print(f'Mean Kuhn length in nm: {np.mean(kuhns_36to47_no_unwrap):.2f}')
#try Andy's formula on this distribution of linker lengths (no unwrapping)
print(np.mean(1/kuhns[35:47])**(-1))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14622980838401675, 'fraction_numerical': 0.055065070956771954, 'mean_word_length': 5.602763385146805, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '119844', 'license': 'MIT', 'n_tokens_mistral': 6774, 'n_tokens_neox': 6037, 'n_words': 1519, 'path': 'notebooks/kinked_WLC_R2_070218.ipynb', 'repo_name': 'brunobeltran/nuc_chain', 'size': '15291'}