Apomyoglobin

---

Summary:

In our previous work (DOI: 10.1002/jcc.23738)

Warning: This notebook will require a significant amount of memory.

---

import numpy as np
import pandas as pd
import biceps
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings('ignore', category=FutureWarning)

BICePs - Bayesian Inference of Conformational Populations, Version 2.0

Warning on use of the timeseries module: If the inherent timescales of the system are long compared to those being analyzed, this statistical inefficiency may be an underestimate.  The estimate presumes the use of many statistically independent samples.  Tests should be performed to assess whether this condition is satisfied.   Be cautious in the interpretation of the data.

####### Data and Output Directories #######
print(f"Possible input data extensions: {biceps.toolbox.list_possible_extensions()}")
T=[0,1,4,9,10,12,14,16,18,19,20,21,24]
states=len(T)
datadir="apomyoglobin/"
top=datadir+'pdb/T1/state0.pdb'
dataFiles = datadir+'new_CS_PF'
input_data = biceps.toolbox.sort_data(dataFiles)
print(f"Input data: {biceps.toolbox.list_extensions(input_data)}")
energies_filename = datadir+'energy_model_1.txt'
energies = np.loadtxt(energies_filename)
energies -= energies.min() # set ground state to zero, just in case
states = len(energies)

Possible input data extensions: ['H', 'Ca', 'N', 'J', 'noe', 'pf']
Input data: ['H', 'Ca', 'N', 'pf']

options = biceps.get_restraint_options(input_data)
pd.DataFrame(options)

	ref	sigma	use_global_ref_sigma	extension	weight	file_fmt	precomputed	pf_prior	Ncs_fi	Nhs_fi	beta_c	beta_h	beta_0	xcs	xhs	bs	states
0	uniform	[0.05, 20.0, 1.02]	True	H	1	pickle	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
1	uniform	[0.05, 20.0, 1.02]	True	Ca	1	pickle	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2	uniform	[0.05, 20.0, 1.02]	True	N	1	pickle	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
3	uniform	[0.05, 20.0, 1.02]	True	pf	1	pickle	False	NaN	NaN	NaN	(0.05, 0.25, 0.01)	(0.0, 5.2, 0.2)	(-10.0, 0.0, 0.2)	(5.0, 8.5, 0.5)	(2.0, 2.7, 0.1)	(15.0, 16.0, 1.0)	NaN

####### Parameters #######
nsteps = 1000000
n_lambdas = 2
outdir = '%s_steps_%s_lam'%(nsteps, n_lambdas)
biceps.toolbox.mkdir(outdir)
lambda_values = np.linspace(0.0, 1.0, n_lambdas)
options = [
        {"ref": 'exponential', "weight": 1/3},
        {"ref": 'exponential', "weight": 1/3},
        {"ref": 'exponential', "weight": 1/3},
        {"ref": 'exponential', "weight": 1, "pf_prior": datadir+'b15.npy',
            "Ncs_fi": datadir+'input/Nc', "Nhs_fi": datadir+'input/Nh', "states": T}
        ]
pd.DataFrame(options)

	ref	weight	pf_prior	Ncs_fi	Nhs_fi	states
0	exponential	0.333333	NaN	NaN	NaN	NaN
1	exponential	0.333333	NaN	NaN	NaN	NaN
2	exponential	0.333333	NaN	NaN	NaN	NaN
3	exponential	1.000000	apomyoglobin/b15.npy	apomyoglobin/input/Nc	apomyoglobin/input/Nh	[0, 1, 4, 9, 10, 12, 14, 16, 18, 19, 20, 21, 24]

####### MCMC Simulations #######
for lam in lambda_values:
    ensemble = biceps.Ensemble(lam, energies)
    ensemble.initialize_restraints(input_data, options)
    sampler = biceps.PosteriorSampler(ensemble)
    sampler.sample(nsteps, verbose=False)
    sampler.traj.process_results(f"{outdir}/traj_lambda{lam}.npz")

100%|█████████████████████████████████████████████████████████| 1000000/1000000 [00:27<00:00, 36768.86it/s]

Accepted 59.564499999999995 %


Accepted [ ...Nuisance paramters..., state] %
Accepted [19.0833 19.0586 19.0551  0.2746  0.2746  0.2746  0.2746  0.2746  0.2746
  0.2746  2.0929] %

100%|█████████████████████████████████████████████████████████| 1000000/1000000 [00:27<00:00, 35870.99it/s]

Accepted 61.9413 %


Accepted [ ...Nuisance paramters..., state] %
Accepted [19.1311 19.0861 19.0675  0.2291  0.2291  0.2291  0.2291  0.2291  0.2291
  0.2291  4.4275] %

####### Posterior Analysis #######
A = biceps.Analysis(outdir, nstates=states)
A.plot(plottype="step")

not all state sampled, these states [1 3 4 9] are not sampled
not all state sampled, these states [ 1  2  3  4  9 12] are not sampled

# NOTE: The following cell is for pretty notebook rendering

from IPython.core.display import HTML
def css_styling():
    styles = open("../../../theme.css", "r").read()
    return HTML(styles)
css_styling()