请你帮我理解一下,在给定的时间内,如何计算单个粒子随机移动的均方位移。我读过很多关于这方面的文章(包括Saxton,1991,单粒子跟踪:扩散系数的分布),但仍然很困惑(没有得到正确的答案)。
让我首先向你们展示我是如何做到的,如果我错了,请纠正我:我的做法如下:
1.运行从t=0到t=100的程序
2.计算每一步的位移(s(t)-s(t+tau))。在t=1,2,3,...100),并将其存储在向量中
3.把答案与数字2打成正方
4.找出答案的平均值
本质上,这就是我在Matlab中所做的
%初始化由16个非零格点组成的方格,然后按以下方式计算MSD:
for t=1:tend
% Allow the particle to move randomly in the lattice. Then do the following
[row,col]=find(lattice>0);
centroid=mean([row col]);
xvec=[xvec centroid(2)];
yvec=[yvec centroid(1)];
k=length(xvec)-1; % Time
dt=1;
diffx = xvec(1:k) - xvec((1+dt):(k+dt));
diffy = yvec(1:k) - yvec((1+dt):(k+dt));
xsquare = diffx.^2;
ysquare = diffy.^2;
MSD=mean(xsquare+ysquare);
end为了计算扩散系数,我试图找到MSD。请注意,我正在建模一组格点(16),以表示单个粒子(更真实的生物),而不是仅仅一个。我对for循环中的评论很简短,因为它很长,但是我很高兴将它发送给您。
到目前为止,我得到的MSD值非常小(在0.001到1之间),而我应该在(10-50)范围内得到值。粒子移动的距离很远,所以我的0.001-1的范围肯定是不对的!
这是一篇文章的摘录,我试图复制他们的形象:
“我们首先在一维中对单个细胞进行了一些模拟。我们允许细胞移动一定数量的蒙特卡罗时间步长(MCS),计算出当时移动的平均平方距离,重复了500次,并对这个t的均方距离进行了评估。然后我们重复了10次,得到了平均值。选择这种重复的原因是为了将运行模拟所需的时间保持在合理的水平上,同时确保平均值的标准差相对较小(<7%)。”
你可以在这里查阅“从离散到连续的生物细胞运动模型,2004年,由Turner等人著,物理评论E”。
任何提示都是非常感谢的。
发布于 2012-03-16 12:58:44
粒子沿多少维移动?
我现在还没有Matlab,但是下面是我如何在一维上这样做的:
% pos is the vector of positions
delta = pos(2:100) - pos(1:99);
meanSquared = mean(delta .* delta);发布于 2012-03-24 13:45:44
首先,为什么有一个粒子覆盖多个晶格点?对于MSD来说,最终重要的是质心的位移,它可以表示为一个点。如果您的粒子(或单元格)很大,或者只是迈出了很大的一步,您总是可以做一个更宽的网格。另外,如果你试图从其他地方复制一个图形,你应该使用同样的算法。
对于你的蒙特卡罗模拟,你是做什么的?如果你真正想要的是得到一个位移,你可以一次生成一串随机运动矢量(使用rand或randi),并使用cumsum计算位置。另外,你是否绘制了随机游走图,以确保数据是合理的?
然后,您的代码看起来有点滑稽(见注释)。为什么不直接使用this answer中提供的代码从位置计算MSD呢?
for t=1:tend
% Allow the particle to move randomly in the lattice. Then do the following
[row,col]=find(lattice>0); %# what do you do this for?
centroid=mean([row col]);
xvec=[xvec centroid(2)];
yvec=[yvec centroid(1)]; %# till here, I have no idea what you want to do
k=length(xvec)-1; % Time %# you should subtract dt here
dt=1; %# dt should depend on t!
diffx = xvec(1:k) - xvec((1+dt):(k+dt));
diffy = yvec(1:k) - yvec((1+dt):(k+dt));
xsquare = diffx.^2;
ysquare = diffy.^2;
MSD=mean(xsquare+ysquare);
endhttps://stackoverflow.com/questions/9737578
复制相似问题