%Plot the prediction for the CF position (xNew) as a function of Bicoid
%dosage (D)
%Start by defining the range of D values we're going to calculate the
%prediction for
%We could do
D=[0.5,1,2];
%However, let's create a set of numbers between reasonable bounds using
%linspace
D=linspace(0.3,5); %Create 100 equally spaced points between 0.3 and 5
%Define the model parameters
xCF=0.32; %Fraction of the embryo at which the WT CF shows up
lambda=0.19; %Bicoid decay length in units of fraction of embryo length
%Calculate xNew
xNew=xCF+lambda*log(D);
%Plot it
semilogx(D,xNew,'-r')
xlabel('Dosage')
ylabel('xNew')
%Now analyze the data so we can test our model
%Load an image of a 1x embryo
Embryo=imread('1xBcd.tif');
imshow(Embryo)
%Strategy for finding the CF: We're going to click along the ventral side
%in order to measure the total length of the embryo and the length from the
%anterior end to the CF.
%We can record clicks on an image using ginput. By calling ginput, you can
%click as many times as you want on the plot. When you're done, you press
%enter. We'll keep track of the click corresponding to the CF.
%Clicks=ginput; %Note that we comment this line out after we run the
%script for the first time in order to avoid having to
%click again
CFClick=10; %The click corresponding to the CF
%See how good of a job we did at tracing the contour of the embryo by
%overlaying the image with a plot of the clicks
%First, extract the x and y values from Clicks
xClicks=Clicks(:,1); %All rows in the first column of Clicks
yClicks=Clicks(:,2); %All rows in the second column of Clicks
%Now, make the plots
hold on
plot(xClicks,yClicks,'.-r')
hold off
%Calculate the length of each segment
%First, create "Delta" vectors that have, for example,
%DeltaX=[x2-x1,x3-x2,...]
%This can be done using the diff function
DeltaX=diff(xClicks);
DeltaY=diff(yClicks);
%Second, calculate the square of each Delta element. Note that we need to
%raise to the power of 2 in an element-wise fashion by using the "."
%operator
DeltaX2=DeltaX.^2;
DeltaY2=DeltaY.^2;
%Third, add DeltaX2 and DeltaY2
SquareLengths=DeltaX2+DeltaY2;
%Finally, take the square root
Lengths=sqrt(SquareLengths);
%Calculate the total axis length by summing over all elements of Lengths
LengthEmbryo=sum(Lengths);
%Calculate the length until the CF by adding Lengths from 1 to CFClick-1
LengthCF=sum(Lengths(1:CFClick-1));
%Get the CF position with respect to the total length of the embryo
LengthCF/LengthEmbryo
%We got the data for the different dosages
DData=[0.5,1,2];
xNewData=[0.25,0.32,0.36];
%Plot the prediction
semilogx(D,xNew,'-r')
xlabel('Dosage')
ylabel('xNew')
hold on
plot(DData,xNewData,'.k','MarkerSize',20)
hold off