MLE ddG errors != calculated ddG errors for open cycle

[jeremyaq]

First of all many thanks for your nice contribution to the field!

I was surprised to see that MLE ddG errors differ from calculated ddG errors for open cycles, this is not expected right? Please find bellow the code to reproduce this unexpected results, thanks in advance for your insights!

# cinnabar 0.4.1
from cinnabar import femap
import numpy as np
import itertools

# dummy calculated and experimental data
experimental_data = {'0': {'dG': -1, 'ddG': 2}, '1': {'dG': -2, 'ddG': 2}, '2': {'dG': -3, 'ddG': 2}, '3': {'dG': -4, 'ddG': 2}}
#calculated_data = {'0->1': {'ligand_i': '0', 'ligand_j': '1', 'dG': 1, 'ddG': 1}, '0->2': {'ligand_i': '0', 'ligand_j': '2', 'dG': 2, 'ddG': 1}, '2->3': {'ligand_i': '2', 'ligand_j': '3', 'dG': 3, 'ddG': 1}}
calculated_data = {'0->1': {'ligand_i': '0', 'ligand_j': '1', 'dG': 1, 'ddG': 2.5}, '0->2': {'ligand_i': '0', 'ligand_j': '2', 'dG': 2, 'ddG': 2.8}, '2->3': {'ligand_i': '2', 'ligand_j': '3', 'dG': 3, 'ddG': 3.1}}


# creating cinnabar csv file
cinnabar_filename = 'cinnabar_input.csv'
with open(cinnabar_filename, 'w') as f:
    f.write("# Experimental block\n")
    f.write("# Ligand, expt_DDG, expt_dDDG\n")
    for entry in experimental_data:
        dG = experimental_data[entry]['dG']
        ddG = experimental_data[entry]['ddG']
        f.write(f"{entry},{dG:.2f},{ddG:.2f}\n")
    f.write('\n')
    f.write('# Calculated block\n')
    f.write('# Ligand1,Ligand2,calc_DDG,calc_dDDG(MBAR),calc_dDDG(additional)\n')
    for entry in calculated_data:
        dG = calculated_data[entry]['dG']
        ddG = calculated_data[entry]['ddG']
        molA = calculated_data[entry]['ligand_i']
        molB = calculated_data[entry]['ligand_j']
        f.write(f"{molA},{molB},{dG:.2f},0,{ddG:.2f}\n")


# create femap
fe = femap.FEMap.from_csv(cinnabar_filename)
# Get MLE generated estimates of the absolute values
fe.generate_absolute_values()
graph = fe.to_legacy_graph()
nodes = graph.nodes(data=True)
node_names = [node[0] for node in nodes]
mle_dG = np.asarray([node[1]["calc_DG"] for node in nodes])
mle_dG_err = np.asarray([node[1]["calc_dDG"] for node in nodes])

# Get calcualted ddGs from MLE-dGs in cinnabar
mle_network = {}
for (a, b), node_pair in zip(itertools.combinations(range(len(mle_dG)), 2), itertools.combinations(node_names, 2)):
    # mle_ddG estimates from MLE-dGs
    mle_ddG = mle_dG[b] - mle_dG[a]
    # mle ddG errors from MLE-dG errors (error propagation)
    mle_ddG_err = (mle_dG_err[b] ** 2 + mle_dG_err[a] ** 2) ** 0.5
    # populate mle_network
    mle_network[f"{node_pair[0]}->{node_pair[1]}"] = {'ligand_i': node_names[a], 'ligand_j': node_names[b], 'dG': mle_ddG, 'ddG': mle_ddG_err}
    #print(node_pair, mle_ddG, mle_ddG_err)

for edge in calculated_data:
    # MLE ddGs and calculated ddGs should be the same as we have an open thermodynamic cycle
    print(f"MLE ddG: {round(mle_network[edge]['dG'], 1)}, Calculated ddG: {round(calculated_data[edge]['dG'], 1)}")
    assert round(calculated_data[edge]['dG'], 1) == round(mle_network[edge]['dG'],1 ), f"MLE ddG and calculated ddG for {edge} are not the same"

    # MLE ddG errors and calculated ddG errors should be the same as we have an open thermodynamic cycle
    print(f"MLE ddG error: {round(mle_network[edge]['ddG'], 1)}, Calculated ddG error: {round(calculated_data[edge]['ddG'], 1)}")
    assert round(calculated_data[edge]['ddG'], 1) == round(mle_network[edge]['ddG'], 1), f"MLE ddG error and calculated ddG error for {edge} are not the same"