%In this script, we will compare theory with experimental data. In class,
%we derived "The French Flag" model for patterning in Drosophila. We'll
%mark the positioning of the cephalic furrow in a WT embryo (Bcd dosage=1)
%and a mutant (Bcd dosage=0.5) by using mouse clicks. Let's begin by
%loading the wildtype embryo image.
WTembryo = imread('wild_type.tif');
imshow(WTembryo, []);
%There are a few things to realize about this image. Firstly, the embryo is
%tilted and the exact position of the cephalic furrow is hard to pin
%point. We will have to use some trigonometry to determine the position of
%the anterior-posterior axis and the relative positioning of the cephalic
%furrow. To begin, we'll click on the position of the anterior, the
%posterior, and the cephalic furrow. Wherever you want to define these
%points is fine so long as you are consistent for each embryo. To record
%the positions of the mouse clicks, we'll use MATLABs 'ginput' command.
%Find the anterior and posterior poles.
WTAPclick = ginput(2); %The "2" tells matlab we'll make two clicks.
WTAx = WTAPclick(1,1); %Extract the x position of the first click.
WTAy = WTAPclick(1,2); %Extract the y position of the first click.
WTPx = WTAPclick(2,1); %Extract the x position of the second click.
WTPy = WTAPclick(2,2); %Extract the y position of the second click.
%Now that we have the locations of the anterior and the posterior poles,
%we'll find the position of the cephalic furrow.
WTFclick = ginput(1); %Tell Matlab we'll make only one click.
WTFx = WTFclick(1); %Extract the x component of the click.
WTFy = WTFclick(2); %Extract teh y component of the click.
%We'll start our analysis by finding the length of the vector between the
%anterior and posterior poles as well as the length between the cephallic
%furrow and the anterior pole.
WTAP = sqrt((WTPy - WTAy)^2 + (WTPx - WTAx)^2); %AP distance.
WTFA = sqrt((WTFy - WTAy)^2 + (WTFx - WTAx)^2); %CFA distance.
%We can use trigonometry to find the length of the fragment of the anterior
%-posterior that intersects with the cephalic furrow. Alpha is an interior
%angle of the triangle with hypotenouse AP. Beta is an interior angle in
%the triangle formed with the hypotenouse CFA. Gamma is an interior angle
%of the triangle with CFA as the hypotenouse and a portion of AP as one
%side.
alpha = atan((WTPy - WTAy) / (WTPx - WTAx));
beta = atan((WTFy - WTAy) / (WTFx - WTAx));
gamma = beta - alpha;
%With this information, we can compute the absolute length (in pixels) of
%the fragment of AP between the anterior and the cephalic furrow as well as
%the relative length.
AbsCF = cos(gamma) * WTFA;
RelCF = AbsCF / WTAP;
%As a sanity check, we know that the cephalic furrow in the WT embryo
%should be at about 30% of the length of the embryo. Let's see what we
%calculated.
RelCF; %Should be close to 0.3.
%So now that we know how to calculate the position, let's do the same
%procedure to the Bicoid mutant. We could simplify this by defining a
%function, but for now let's rewrite the same protocol.
%Read the bicoid mutant embryo
BcdEmbryo = imread('bcd_mut.tif');
imshow(BcdEmbryo, []);
Bcdclick = ginput(2);
BcdAx = Bcdclick(1,1);
BcdAy = Bcdclick(1,2);
BcdPx = Bcdclick(2,1);
BcdPy = Bcdclick(2,2);
%Find the cephalic furrow
BcdFclick = ginput(1);
BcdFx = BcdFclick(1);
BcdFy = BcdFclick(2);
%Compute angles and lengths.
alpha = atan((BcdPy - BcdAy) / (BcdPx - BcdAx));
beta = atan((BcdFy - BcdAy) / (BcdFx - BcdAx));
gamma = beta - alpha;
%Anterior-posterior axis.
BcdAP = sqrt((BcdPy - BcdAy)^2 + (BcdPx - BcdAx)^2);
BcdFa = sqrt((BcdFy - BcdAy)^2 + (BcdFx - BcdAx)^2);
%Absolute CF length.
BcdAbsCF = cos(gamma) * BcdFa;
BcdRelCF = BcdAbsCF / BcdAP;
%Now we have the positioning of the cephalic furrow of the WT and the
%bicoid mutant. Let's compare this data with the theory. We'll compute the
%curve we defined in class given our WT cephalic furrow position and plot
%our data points on the same set of axes. We'll start by defining some
%parameters.
lambda =0.2; %Decay length of Bicoid (unitless).
xcf = RelCF; %Position of the furrow on a standard embryo of length 500µm.
Dosages = linspace(0.5, 3); %Range of Bcd dosages.
NewXcf = xcf + log(Dosages)*lambda; %Our theoretical curve.
%Now plot it!
semilogx(Dosages, NewXcf);
hold on
plot(1, xcf, 'ro'); %Plot the WT as open red circles.
plot(0.5,BcdRelCF, 'ko'); %Plot the WT as open black circles.
xlabel('Bicoid dosage (fractional)')
ylabel('Cephalic furrow position from anterior (relative units)');
legend('Prediction', 'WT Position', 'Mutant position', 'Location', 'Northwest');
hold off
%How does the data compare to the model? Does the data have to match the
%theory?