Understanding the output

Running a ligandswap calculation is (hopefully) straightforward. Analysis is more involved…

The output of the ligandswap calculation is written into the directory output. This directory contains summaries of the analysis performed after each iteration of sampling, and PDB files that contain the coordinates of the atoms.

Of most interest is the free energy that is actually calculated by ligandswap. All of the free energy statistics that are needed to calculate the free energy are stored in the file output/freenrgs.s3. To get this, we need to use the analyse_freenrg program that comes with Sire. This program is used to analyse the free energy output from many of the Sire-based applications.

To use analyse_freenrg, please type;

$SIRE/bin/analyse_freenrg -i output/freenrgs.s3

The -i option tells analyse_freenrg the name of the s3 file that contains the free energy components to analyse.

The output of this analysis is printed to the screen.

The first part of the output is shown below;

# Convergence
# Iteration 
# Bennetts # FEP # TI 
1 518.0732796078831 518.7156698068001 519.6264675640966 
2 513.3134777256312 515.4536030784152 514.1634743397709 
3 513.9175587894974 515.6958276498924 514.812530912826 
4 513.2527480088207 515.4313733065027 514.0776866026532 
5 511.47566149724435 513.5609397154682 511.28529099672755 

Plot of free energy versus iteration
   │                                                                            
   ┼+513.911      .             .              .             .                  
   │                            .                            .              .   
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   ┼+10.9122                                                                    
───┼┼───────────────────────────────────────────────────────────────────────┼───
   │ +0.026                                                         +4.98688    

This shows the relative binding free energy calculated using the statistics collected only during each iteration of Monte Carlo sampling. In our case, we performed five iterations. This means that we have five separate sets of estimates of the relative binding free energy. ligandswap calculates the relative binding free energy using the Bennett Acceptance Ratio, free energy perturbation (FEP) and thermodynamic intergration (TI) methods, meaning that each iteration results in three independent estimates of the relative binding free energy.

Below this, is a graph of the free energy predicted at each iteration, versus iteration, which we call the “convergence graph”. As we have only run a (very) short calculation, the estimate is high, and there is little variation. As we will see later, for longer simulations, this convergence graph can give us good indication as to whether or not the free energy prediction has converged.

Next, analyse_freenrg outputs the potentials of mean force (PMFs) across λ predicted by each of the three free energy methods, e.g. here is the Bennetts (BAR) prediction;

# Bennetts
# Lambda  PMF  Maximum  Minimum 
0.0  0.0  0.0  0.0
0.005  -14.140560913267471  -14.13819614006264  -14.142925686472303
0.071  -154.97080622520267  -154.86387728229295  -155.07773516811238
0.137  -237.4119541405933  -237.027032095831  -237.79687618535561
0.203  -287.4699709034169  -287.0745659012938  -287.86537590553996
0.269  -317.5955934182356  -317.1065931464901  -318.08459368998103
0.335  -334.67538292822906  -334.15764241188975  -335.19312344456836
0.401  -342.5482903205475  -342.02729130610203  -343.06928933499296
0.467  -343.37768089490476  -342.74426418577633  -344.0110976040332
0.533  -337.64691807126593  -336.92740610595376  -338.3664300365781
0.599  -323.4608847456288  -322.50692020015595  -324.41484929110163
0.665  -296.35848200242947  -295.4001026136855  -297.31686139117346
0.731  -251.80059916468298  -250.8415289453025  -252.75966938406344
0.797  -185.44235518658832  -184.38297304774196  -186.50173732543468
0.863  -76.09708201207764  -74.77671322745805  -77.41745079669722
0.929  122.9104621765241  124.23269864588933  121.58822570715888
0.995  469.12662131515594  470.4488577845212  467.8043848457907
1.0  504.5267693078106  505.86002499289117  503.19351362273

PMF Plot of free energy versus lambda
   │                                                                            
   ┼+495.2                                                                  ○   
   │                                                                       ○    
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                   ○        
   │                                                                            
───○┼───────────────────────────────────────────────────────────────────────┼───
   │ +0.0051                                                              +1    
   │                                                              ○             
   │    ○                                                                       
   │                                                         ○                  
   │         ○                                          ○                       
   ┼-325.572      ○    ○                       ○   ○                            
   │                       ○    ○    ○    ○                                     
   │                                                                            

here is the FEP prediction;

# FEP
# Lambda  PMF  Maximum  Minimum 
0.0  0.0  0.0  0.0
0.005  -14.140560913267471  -14.13819614006264  -14.142925686472303
0.071  -155.0031405345896  -154.79692672764497  -155.2093543415342
0.137  -237.89239377740927  -236.57335494516124  -239.2114326096573
0.203  -287.87714680838593  -286.3002680211615  -289.45402559561035
0.269  -317.94050073450995  -316.1684209898753  -319.7125804791446
0.335  -335.06177709734925  -333.13817573571436  -336.98537845898414
0.401  -342.95029089079213  -340.98139650775136  -344.9191852738329
0.467  -343.61872153009654  -341.0662129283602  -346.1712301318329
0.533  -337.4047187205829  -333.58748310706557  -341.2219543341002
0.599  -322.59982677056627  -317.31157116967734  -327.8880823714552
0.665  -295.50238754071313  -290.18828737369563  -300.81648770773063
0.731  -251.35043029735021  -244.7850864513008  -257.91577414339963
0.797  -182.09274165262335  -168.30119518833135  -195.88428811691534
0.863  -68.14262923130357  -43.98837648596146  -92.2968819766457
0.929  130.8743097875974  155.05749744042296  106.69112213477187
0.995  477.1185724821613  501.37254734192595  452.8645976223967
1.0  512.518720474816  536.7837145502959  488.25372639933596

PMF Plot of free energy versus lambda
   │                                                                            
   ┼+503.101                                                                ○   
   │                                                                       ○    
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                            
   │                                                                   ○        
   │                                                                            
───○┼───────────────────────────────────────────────────────────────────────┼───
   │ +0.0051                                                              +1    
   │                                                              ○             
   │    ○                                                                       
   │                                                         ○                  
   │         ○                                          ○                       
   ┼-325.64       ○    ○                       ○   ○                            
   │                       ○    ○    ○    ○                                     
   │                                                                            

and here is (an extract of) the TI prediction;

# TI
# Lambda  PMF  Maximum  Minimum 
0.0  0.0  0.0  0.0
0.0  -27.690011993014004  -26.87010722066145  -28.509916765366558
0.01  -53.316860078252446  -51.76275284604756  -54.87096731045733
0.02  -76.97371007870679  -74.83334429070639  -79.11407586670718
0.03  -98.77365387616248  -96.21397662399913  -101.33333112832584
0.04  -118.84123556325576  -116.01982891640199  -121.66264221010954
====
0.9700000000000006  381.30811192051385  398.6528636778263  363.9633601632014
0.9800000000000006  444.33480229691247  462.3265817367379  426.343022857087
0.9900000000000007  513.5784295374591  530.8436008461455  496.3132582287727
1.0000000000000007  513.5784295374542  530.8436008461407  496.31325822876767

PMF Plot of free energy versus lambda
   │                                                                            
   │                                                                       ○○   
   ┼+504.133                                                                    
   │                                                                      ○     
   │                                                                      ○     
   │                                                                     ○      
   │                                                                    ○       
   │                                                                   ○        
   │                                                                   ○        
   │                                                                 ○○         
   │                                                                ○           
───○┼──────────────────────────────────────────────────────────────○○───────┼───
   ○○+0.0051                                                      ○       +1    
   │ ○○                                                        ○○○              
   │   ○○                                                    ○○○                
   │     ○○○                                              ○○○                   
   │        ○○○○                                      ○○○○○                     
   │            ○○○○○○○○                     ○○○○○○○○○                          
   ┼-327.065            ○○○○○○○○○○○○○○○○○○○○○                                   
   │                           

As you can see, a graph of the PMF (free energy versus λ) is printed below the actual numerical values. The graph is useful to show you the shape of the PMF. The numerical values also give estimates of the error, i.e. for TI, the relative free energy is predicted to be between 496 kcal mol-1 and 531 kcal mol-1 (all energies are in kilocalories per mole). These error ranges are always massive overestimates of the statistical error for FEP and TI, and massive underestimates for BAR.

In my opinion, the best way to get an estimate of the error is to look at the level of agreement between the different methods of calculating the relative free energy. This is summarised at the end of the output of analyse_freenrg and is shown below;

# Free energies 
# Bennetts = 504.5267693078106 +/- 1.3332556850805601 kcal mol-1#
# FEP = 512.518720474816 +/- 24.264994075479983 kcal mol-1#
# TI = 513.5784295374542 +/- 17.26517130868654 kcal mol-1 (quadrature = 528.9284057323131 kcal mol-1)#

This shows four estimates of the relative binding free energy, together with their errors. There is about a 25 kcal mol-1 disagreement between the methods, suggesting that the calculation has not really converged. Converged calculations should typically agree to within 0.5 to 1.0 kcal mol-1.

Now, you may ask, why do we have four free energy estimates from three methods? This is because there are two ways in which we can obtain the relative binding free energy from TI; quadrature using the trapezium rule, and analytic integration of a polynomial that is automatically fitted to the gradients. The main TI value (and the curve that is plotted) show the result of analytic integration, while the quadrature result is printed as “quadrature”. analyse_freenrg prints both, as the quadrature result is much more sensitive to poor convergence, and so poor agreement of this with the other methods is a really good indication of inherent error.


Previous Up Next