3 Use Voronoi polygons to compute the lifetime of lipids near the surface of
7 Alan Grossfield, University of Rochester Medical Center
20 Computes the autocorrelation function for a set of timeseries.
22 Destroys the timeseries matrix.
23 Does not trap for timeseries with zero std dev (e.g. all constant)
ave = numpy.mean(timeseries, axis=1, keepdims=True
dev = numpy.std(timeseries, axis=1, keepdims=True
length = timeseries.shape
corr = numpy.zeros([expected, len(timeseries)], float)
num = timeseries[:,i:] * timeseries[:,:num_vals]
corr[i] = numpy.average(num, axis = 1)
ave_corr = numpy.average(corr, axis=1)
dev_corr = numpy.std(corr, axis=1)
49 Usage: lipid_lifetime.py system trajectory zmin zmax protein-selection all-lipid-selection target-lipids
51 This program uses 2D Voronoi decomposition to identify lipids in contact
52 with the protein surface, and returns the time-correlation function
53 for this contact (averaged over all lipids considered). This can be fit
54 to a sum of exponentials to extract a "lifetime".
56 zmin and zmax define the "slice" of the membrane and protein used in the
57 Voronoi decomposition. The details of the choice shouldn't matter too much,
58 as long as you pick out one leaflet or the other; beyond that, your choice
59 should be guided by the shape of the protein and the variation of whatever atom
60 you use to define the lipid. At the moment, a lipid that leaves the z-slice
61 is treated the same as one that leaves the protein, so I'd say handling of
62 flip-flop is not very robust.
64 protein-selection is the set of atoms used to define the "protein", or the
65 thing that you're calculating the lifetime around. It is assumed to be 1 thing
66 (e.g. all the atoms are treated together, and we actually just use the
67 centroid), so you couldn't put a selection for something like a bunch of
68 cholesterol molecules there and expect to get the lifetime for lipids around
69 cholesterol. A fairly sparse selection is probably best, e.g. just alpha
72 all-lipids is a selection that described all of the lipids that are found in
73 the leaflet you're examining. The code is written assuming that you'll pick
74 one atom per lipid, e.g. the phosphorus. The important thing is that all
75 lipids (and cholesterols, etc) in the leaflet in question must be represented,
76 or else the resulting "gaps" may cause some lipids to be incorrectly be
77 marked as neigboring the protein.
79 target-lipids is a selection of those lipids for which you want to compute the
80 lifetime. This group is by definition a subset of the all-lipids selection;
81 since we actually apply this selection to the group created by the all-lipids
82 selection, you can use 'all' as this selection if you want to compute the
83 correlation function for every lipid. Alternatively, if you have a mixture
84 of POPE and POPC lipids, but you want just the lifetime for just the POPE
85 lipids, you might use 'resname =~ "POP[CE]" and name == "P"' for all-lipids
86 and 'resname == "POPE"' for target-lipids.
89 Note: The third column in the output is the standard deviation in the
90 correlation function, where you've averaged all lipids. If you believe
91 that the lipids can be treated as independent measurements, you would divide
92 this value by sqrt(N_lipids_considered) to get the standard error in the
93 correlation function, which is related to the statistical uncertainty. I
94 didn't do this because it's not clear to me that that is a good assumption.
__name__ == '__main__'
len(sys.argv) <= 1 or
sys.argv == '--fullhelp'
103 print "#"
, " "
system_file = sys.argv
traj_file = sys.argv
zmin = float(sys.argv)
zmax = float(sys.argv)
protein_selection = sys.argv
all_lipid_selection = sys.argv
target_lipid_selection = sys.argv
slicer = ZSliceSelector(zmin, zmax)
neighbor_timeseries = numpy.zeros([len(target_lipids), len(traj)], numpy.float)
protein_slice = slicer(protein)
centroid = protein_slice.centroid()
all_lipids_slice = slicer(all_lipids)
target_lipid_slice = slicer(target_lipids)
combined = protein_centroid + all_lipids_slice
v = VoronoiWrapper(combined, padding)
protein_region = v.get_region_from_atomid(protein_centroid.id())
region = v.get_region_from_atomid(target_lipid_slice[i].id())
neighbor_timeseries[i][traj.index()] = 1
ave = numpy.mean(neighbor_timeseries, axis=1)
is_nonzero = numpy.select([abs(ave)>1e-6], [ave])
neighbor_timeseries = numpy.compress(is_nonzero,
ave_corr, dev_corr = autocorrel(neighbor_timeseries, neighbor_timeseries.shape/2)
179 print "#Frame\tCorrel\tDev"
i, ave_corr[i], dev_corr[i]
Python-based wrapper for LOOS Trajectories This class turns a loos Trajectory into something more pyt...
Basic Atom class for handling atom properties.
AtomicGroup selectAtoms(const AtomicGroup &source, const std::string selection)
Applies a string-based selection to an atomic group...
Class for handling groups of Atoms (pAtoms, actually)
AtomicGroup createSystem(const std::string &filename, const std::string &filetype)