### #77 Brownian Carnot engine-I 微觀卡諾引擎-I

Hello. Today we will use MATLAB program for simulation of gradient sensing by slime mold. In our previous episode we have write down explicitly the estimated concentration gradient in z direction by monitoring the position of inward particle flow and summing their cosine. That is

So let's see how efficient and accurate it is.

## Our MATLAB script

%% Basic setting
tic; clear; clc;
r_Dicty = 0;
location_Dicty = [10,0, 0];
x_source = 0;
width_source = 15;
x_Boundary = 20;
speed = 0.1;
D = speed^2/6;
t = 1;
tRun = 40000;
Detected_List = nan(tRun, 3);
x_particle_all = nan(tRun, 3);
v_particle_all = nan(tRun, 3);
location_Dicty_rpt = repmat(location_Dicty, tRun, 1);


In the first part of our code we simply set the location of our slime mold cell and the speed of our nutrient particle. The diffusion constant D is calculated by knowing in 3D world, <r^2>=6DT.

%% get into steady state before our experiment
while t<=tRun
theta = rand(1,1)*pi; phi = rand(1,1)*2*pi;
y_source = width_source*(rand(1,1)-0.5);
z_source = width_source*(rand(1,1)-0.5);
x_particle_all(t, :) = [x_source, y_source, z_source];
v_particle_all(t, :) = [cos(phi)*sin(theta), sin(phi)*sin(theta), cos(theta)];
%move
x_particle_all = x_particle_all + v_particle_all*speed;
%check collapse boundary
a = x_particle_all(:,1)<0 -="" a="" b="x_particle_all(:,1)" v_particle_all="" x_particle_all=""> x_Boundary;
x_particle_all(b,:) = nan;
v_particle_all(b,:) = nan;
c = x_particle_all(:,2) > 0.5*width_source;
x_particle_all(c,2) = x_particle_all(c,2) - width_source;
d = x_particle_all(:,2) < -0.5*width_source;
x_particle_all(d,2) = x_particle_all(d,2) + width_source;
e = x_particle_all(:,3) > 0.5*width_source;
x_particle_all(e,3) = x_particle_all(e,3) - width_source;
f = x_particle_all(:,3) < -0.5*width_source;
x_particle_all(f,3) = x_particle_all(f,3) + width_source;
%check collapse detector
dist_all = (sum((x_particle_all - location_Dicty_rpt).^2, 2)).^0.5;
e = dist_all <= r_Dicty;
Detected_List(e,:) = x_particle_all(e,:);
x_particle_all(e,:) = nan; v_particle_all(e,:) = nan;
%change direction
b = ~isnan(v_particle_all);
c = b(:,1)&b(:,2);
d = 2*pi*rand(tRun,1); e = pi*rand(tRun,1);
v_particle_all = [sin(e).*cos(d), sin(e).*sin(d), cos(e)];
v_particle_all(~c,:) = nan;
%next round
t = t + 1;
%display(t);
end


Because our experimental environment is a 15x15x20 hexahedral, which is quite large and is able to accommodate lots of particles, so before our experiment we must have our particle diffusion get into steady state. The above code could also be written in object-oriented coding. However, it would be much much more efficient in MATLAB if you could express everything you need to do in matrix form. The particles are generated in the left boundary and disappear if they collapse the right boundary. Let's see the result of above code:

Hmm, it looks like it is now in steady state. (Of course there are some ways to check if it is in steady state already. How will you do?) So now let's place our slime mold in the center.

%% place our slime mold cell
r_Dicty = 3;
dist_all = (sum((x_particle_all - location_Dicty_rpt).^2, 2)).^0.5;
e = dist_all <= r_Dicty;
x_particle_all(e,:) = nan; v_particle_all(e,:) = nan;
x_particle_all(isnan(x_particle_all(:,1)),:)=[];
v_particle_all(isnan(v_particle_all(:,1)),:)=[];
nSurvived = length(x_particle_all(:,1));
t = nSurvived + 1;
tRun = 60000;
location_Dicty_rpt = repmat(location_Dicty, tRun, 1);
x_particle_all_2 = nan(tRun, 3);x_particle_all_2(1:nSurvived, :) = x_particle_all;
v_particle_all_2 = nan(tRun, 3);v_particle_all_2(1:nSurvived, :) = v_particle_all;
x_particle_all = x_particle_all_2;
v_particle_all = v_particle_all_2;


And we could start our experiment:

%% start experiment
tMeasure = 300;
nMeasure = ((tRun - nSurvived) - mod((tRun - nSurvived), tMeasure))/tMeasure;
cz_estimated = zeros(nMeasure,1);
while t<=tRun
theta = rand(1,1)*pi; phi = rand(1,1)*2*pi;
y_source = width_source*(rand(1,1)-0.5);
z_source = width_source*(rand(1,1)-0.5);
x_particle_all(t, :) = [x_source, y_source, z_source];
v_particle_all(t, :) = [cos(phi)*sin(theta), sin(phi)*sin(theta), cos(theta)];
%move
x_particle_all = x_particle_all + v_particle_all*speed;
%check collapse boundary
a = x_particle_all(:,1)<0 -="" a="" b="x_particle_all(:,1)" v_particle_all="" x_particle_all=""> x_Boundary;
x_particle_all(b,:) = nan;
v_particle_all(b,:) = nan;
c = x_particle_all(:,2) > 0.5*width_source;
x_particle_all(c,2) = x_particle_all(c,2) - width_source;
d = x_particle_all(:,2) < -0.5*width_source;
x_particle_all(d,2) = x_particle_all(d,2) + width_source;
e = x_particle_all(:,3) > 0.5*width_source;
x_particle_all(e,3) = x_particle_all(e,3) - width_source;
f = x_particle_all(:,3) < -0.5*width_source;
x_particle_all(f,3) = x_particle_all(f,3) + width_source;
%check collapse detector
dist_all = (sum((x_particle_all - location_Dicty_rpt).^2, 2)).^0.5;
e = dist_all <= r_Dicty;
Detected_List(e,:) = x_particle_all(e,:);
x_particle_all(e,:) = nan; v_particle_all(e,:) = nan;
%change direction
b = ~isnan(v_particle_all);
c = b(:,1)&b(:,2);
d = 2*pi*rand(tRun,1); e = pi*rand(tRun,1);
v_particle_all = [sin(e).*cos(d), sin(e).*sin(d), cos(e)];
v_particle_all(~c,:) = nan;
%next round
if mod((t - nSurvived), tMeasure) == 0
j = (t - nSurvived)/tMeasure;
Detected_List_temp = Detected_List;
Detected_List_temp(isnan(Detected_List_temp(:,1)),:)=[];
nDetected = length(Detected_List_temp(:,1));
location_Dicty_rpt_temp = repmat(location_Dicty,nDetected, 1);
diff_vec_temp = Detected_List_temp - location_Dicty_rpt_temp;
cosine_detect = diff_vec_temp(:,1)./((diff_vec_temp(:,1).^2+diff_vec_temp(:,2).^2+diff_vec_temp(:,3).^2).^0.5);
cz_estimated(j) = sum(cosine_detect)/(4*pi*D*r_Dicty^2*tMeasure);
display(cz_estimated(j));
Detected_List = nan(tRun, 3);
end
t = t + 1;
end


So now our cute slime mold is measuring the gradient for us. However, before we ask him to print the result, we must have some idea about what we expect to get. Assume the ΔT means 1sec for each round. The experimental setting could be formulated into the following problem:

We have simplified our equation based on symmetry of our system already. The solution is very simple. And the estimated concentration gradient in x>0 would be:

So let's see if our expectation meets what we measured:

cz_theoretical = -(1/(width_source*width_source))/(2*D);
cz_best_estimate = mean(cz_estimated);
var_cz_theoretical = -cz_theoretical*x_Boundary/(12*pi*D*r_Dicty^3*tMeasure);
var_cz_estimate = var(cz_estimated);
display(cz_theoretical);
display(cz_best_estimate);
display(var_cz_theoretical);
display(var_cz_estimate);
toc;


And the result we get:

The estimation of variance could be found in suggested reading. They perfectly match each other. We could also display our result in figure:

%% plot result
[x,y,z] = sphere;
figure
surf(r_Dicty*x+location_Dicty(1),r_Dicty*y+location_Dicty(2),r_Dicty*z+location_Dicty(3));
hold on;
axis square;
Detected_List(isnan(Detected_List(:,1)),:)=[];
x_particle_all(isnan(x_particle_all(:,1)),:)=[];
if ~isempty(Detected_List)
scatter3(Detected_List(:,1),Detected_List(:,2),Detected_List(:,3) , 1);
end
scatter3(x_particle_all(:,1),x_particle_all(:,2),x_particle_all(:,3) , 1);
toc;


The tiny blue dots are the recorded particles. Noted that our slime mold only record a little particles compared to the number of particles in the background. It's really amazing, isn't it? And that's all about our this series, gradient sensing. I hope you enjoy it.

/*---------Divider---------*/