%Dependence Shaded area error bar plot %(https://www.mathworks.com/matlabcentral/fileexchange/58262-shaded-area-error-bar-plot) titlename='dhead:dabdomen:dwing=1200:600:400'; %%%%%%%%%%%%% Set up the parameters %%%%%%%%%%%%%%% fly_num=10; % Run the simulation for 10 times d_body=[1200;600;400]; % initial dust on body parts (head,abdomen,wing) d_leg=200; % initial dust on the legs dr_dustremoval=0.004; % how many percent of dust is transfered to legs from a body part in each cleaning iteration %for each leg rubbing, 10*dr_dustremoval% dust is removed ethogram_time_sum=[]; %ethogram for 10 simulation dustamount_time_sum=[]; for i=1:fly_num [dustamount_time_sum(:,:,i),ethogram_time_sum(i,:)]=grooming_model(d_body,d_leg,dr_dustremoval); end %%%%%%%%%%%%% Anterior grooming probability over time %%%%%%%%%%%%%%% AP_trend_window=1000; AP_trend_step=fix(2.*length(ethogram_time_sum)/AP_trend_window-1); A_trend=zeros(2,AP_trend_step); for i=1:size(ethogram_time_sum,1) for j=1:AP_trend_step A_trend(i,j)=(sum(ethogram_time_sum(i,1+(AP_trend_window/2)*(j-1):(AP_trend_window/2)*(j+1))==1)+sum(ethogram_time_sum(i,1+(AP_trend_window/2)*(j-1):(AP_trend_window/2)*(j+1))==2))/AP_trend_window; end end %%%%%%%%%%% Grooming progression figure %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure; x=0:(AP_trend_step-1); plot_areaerrorbar(A_trend,'r','sem',1.8,'--'); hold on; plot_areaerrorbar(1-A_trend,'b','sem',1.8,'--'); hold off; legendlines = findobj(gca,'Type','line'); legend([legendlines(2) legendlines(1)],'Anterior','Posterior'); title(titlename,'FontSize',16); xlabel ('Minutes','FontSize',14); ylabel ('behavior (%active time)','FontSize',14); set(gca,'XTick',[0:7.1:99.4]); set(gca,'XTickLabel',{'0','2','4','6','8','10','12','14','16','18','20','22','24','26','28'}); set(gca,'YTick',[0:0.1:1]); set(gca,'YTickLabel',{'0','10','20','30','40','50','60','70','80','90','100'}); ylim([0 1]); xlim([0 50]); %%%%%%%%%%% Ethogram figure %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure; imagesc(ethogram_time_sum); title(titlename,'FontSize',16); xlabel ('Frames','FontSize',14); ylabel ('Flies','FontSize',14); color = [0.9, 0.6, 0.3;0.4,0,0.9;0, 0.4, 0.0;0, 0.9, 0.0;0, 0.6, 1.0;]; colormap(color); c=colorbar('Ticks',[1.4,2.2,3,3.8,4.6],'TickLabels',{'Front leg','Head','Abdomen','Back Leg','Wing'},'Direction','reverse'); %%%%%%%%%%% Main function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [dustamount_time,etho]=grooming_model(body_dust,leg_dust,clean_incr1) dr=clean_incr1; iteration=2000; %grooming iteration number d_record=zeros(5,iteration); a_record=zeros(5,iteration); etho_record=zeros(1,iteration); d_record(:,1)=[body_dust;leg_dust;leg_dust]; for i=1:iteration for j=1:5 %activity level for head,abdomen,wing,front legs,back legs a_record(j,i)=abs(normrnd(d_record(j,i),d_record(j,i)/5)); end %for winner, 1 means head, 2 means abdomen, 3 means wing, 4 means %front legs, 5 means backlegs [a_level,winner] = max(a_record(:,i)); if i==iteration break; end d_record(:,i+1)=d_record(:,i); d_delta=dr*d_record(winner,i); if winner==4 || winner==5 d_record(winner,i+1)=d_record(winner,i+1)-10*d_delta; else d_record(winner,i+1)=d_record(winner,i+1)-d_delta; end %if head wins, transfer dust to front legs if winner==1 d_record(4,i+1)=d_record(4,i+1)+d_delta; end %if abdomen or wing wins, transfer dust to back legs if winner==2 || winner==3 d_record(5,i+1)=d_record(5,i+1)+d_delta; end %for etho_record, 1 means front legs, 2 means head, 3 means %abdomen, 4 means back legs, 5 means wing (same as normal ethograms) switch winner case 1 etho_record(i)=2; case 2 etho_record(i)=3; case 3 etho_record(i)=5; case 4 etho_record(i)=1; case 5 etho_record(i)=4; end end %%%%%%%%%%%%% Convert grooming events into ethograms (frames)%%%%%%%%%%%%%%% %the distribution of grooming bout duration is from human labeling data %from two dusted CantonS frame_num=25000; etho=zeros(1,frame_num); dustamount_time=zeros(5,frame_num); k=1; for j=1:length(etho_record) switch etho_record(j) case 1 boutlength=gamrnd(1.74571,32.9363); case 2 boutlength=gamrnd(2.80388,8.36872); case 3 boutlength=gamrnd(2.20865,9.68763); case 4 boutlength=gamrnd(3.26172,10.5877); case 5 boutlength=gamrnd(2.17277,8.4282); end boutlength=uint32(boutlength); if k+boutlength>=frame_num etho(k:frame_num)=etho_record(j); dustamount_time(:,k:frame_num)=repmat(d_record(:,j),1,frame_num-k+1); break; else etho(k:k+boutlength-1)=etho_record(j); dustamount_time(:,k:k+boutlength-1)=repmat(d_record(:,j),1,boutlength); k=k+boutlength; end end end