Help data interpretation / steady state

[Victor Montal Blancafort]

Dear Westpa experts, 

I have been running some WE simulations to capture the flux of protein-protein dimerization process. Following your published tutorial, before going forward for fancy analyses, I wanted to ensure that my simulation has converged. This is where my doubts begin. 

I have three states: bound, pre-bound and unbound defined at my .cfg.  (attached files). After running 300+ iterations (4 workers per bin) I stopped my simulation to explore if it had already converged or not.

Using w_ipa -oa, I have generated my direct.h5 file. When I am visualizing the "target_flux_evolution", I am observing a weird spike at late iterations (flux_evolution.png), whereas i was expecting a smooth almost-flat line...  I have also plotted the arrivals values (arrivals.png), and the pattern is not that different along the simulation, hence I suspect the second-spike is the result of a "high-weight" walker reaching the bound state (run trace on them and it seems to trace-back to iteration 1)? 

How should I proceed? More iterations to try to reach a plateu on my flux_evolution? Should I move to haMSM to facilitate convergence?

Thanks in advance for any insight!

KR, 

Victor Montal

[Jeremy Leung]

Hi Victor,

Welcome! (Sorry for the late reply, for some reason I didn't get an email notification for this)

To be fair, your y-axis ticks for the target flux evolution isn't that wide and some fluctuations are not uncommon (I would assume it's reasonably flat at log-scale). But just to be sure, overall I would recommend running a little longer (~50-100 iterations) but there is something you can plot first. The arrivals plot only show the number of segments going in, so it does not completely rule out the larger-weight walkers like you hoped, btw.

Currently, you're looking at target flux with three states defined. And while you do have recycling enabled, I do not know what your WE tau is and how often you're saving your frames (tau/25 frames= ???). The target flux dataset catches segments that jump into the target state, back out & back in all within a single tau. This could lead to double counting and thus a small jump in flux. You would want to recalculate only with two states and use the `conditional_flux_evolution` dataset from direct.h5, and normalize it with the cumulative average of `labeled_population` from assign.h5.

e.g., For state 0 --> state 1:

```

numerator = direct.h5['conditional_flux_evolution'][:,0,1]['expected']  # already cumulatively averaged by w_ipa/w_direct

temp = numpy.sum(assign.h5['labeled_populations'][:, 0], axis=1)  # total probability of segments who last visited 0, summing all bins in each iteration

denominator = numpy.cumsum(temp) / np.arange(1, len(temp)+1)  # cumulative average of the labeled population.

flux = numerator / denominator   # This is in terms of probability, so you'll need to divide by WE tau (e.g. 100 ps = 1e-10 s) to get it in terms of actual molecular time

```

Might be helpful:

https://github.com/westpa/westpa/wiki/HDF5-File-Organization-of-Simulation-Data#user-content-Overall_structure_of_directh5

https://github.com/westpa/westpa/wiki/Calculating-Rate-Constants

Hope this helps,

Jeremy L.