program pullTitin;   {compile in TurboPascal and run on a PC with Turbo
Pascal units (e.g. Graph.tpu) and Borland graphic interface files
(e.g. egavga.bgi) in the default directory}
{$N+}  {uses numeric coprocessor in some PCs}
uses crt,graph;
const
kT = 4.14;    {thermal energy in pN*nm}
A=2.0;        {polypeptide persistence length in nm}
dt=0.02;      {time between iterations in seconds}
speed=65;     {pipette speed in nm/sec.}
trapStif=0.1; {trap stiffness in pN/nm}
fMax=40;      {max force we can pull in pN}
Nmax=70;      {maximum number of native domains}
L=5000;       {maximum contour length in nm}
natLength=4;  {native domain length in nm}
unfLength=35; {unfolded domain length in nm}

var i,j,k,xMax,yMax,direction,N,oldN:integer;
    F,ext,pipetteX,trapX,trapAdjustment,L1,L2,dN:real;
    endOfCycle:boolean;
    ch:char;  message:string;

procedure initializeScreen;
begin
i:=detect;
initgraph(i,j,'');
xMax:=getmaxX; yMax:=getMaxY;
setColor(white);
rectangle(0,0,xMax,yMax);
for i:= 1 to 10 do  {mark x-axis every micron}
   begin j:=round(xMax*i*(1000/L));
   line(j,yMax-10,j,yMax); line(j,0,j,10);
   str(i,message);outTextXY(j-1,20,message+'um');end;
for i:= 1 to fMax div 10 do  {mark y-axis every 10 pN}
   begin j:=yMax-round(yMax*i*(10/fMax));
   for k:=2 to 30 do putPixel(k*xMax div 30,j,white);
   str(i*10,message); outTextXY(10,j-3,message+'pN');
   end;
outTextXY(xMax-250,yMax-20,'q-quit    r-repeat  e-erase');
end;

procedure movePipette;
var dx:real;
    oldDirection:integer;
begin
  oldDirection:=direction;
  if (F>fMax) or (ext>L) then direction:= -1;
  if (F<0) or (ext<L1+10)    then direction := 1;
  dx:=speed*dt;
  pipetteX := pipetteX + direction * dx;
  ext:=pipetteX-trapX;
  endOfCycle:=(direction=1) and (oldDirection=-1);
end;

procedure adjustForce;
var x:real;
begin
x:=(ext-L1)/L2;                  {fractional extension of the unfolded region}
F:=kT*(x+1/(4*sqr(1-x))-0.25)/A;  {Marko/Siggia force function for WLC}
end;

procedure adjustTrap;
var dF,stiffness:real;
begin
if F<0 then F:=0;
dF:=F-trapX*trapStif;                      {force imbalance on trapped bead}
stiffness:=trapStif+(4*F*sqrt(F*A/kT)/L2); {combined trap/titin stiffness}
trapAdjustment:=dF/stiffness;
trapX:=trapX+trapAdjustment;
ext:=pipetteX-trapX;
end;

procedure unfold; {reduces number of folded domains according to force, time}
const Afreq=1E8;  {attempts per second}
      dXunf=0.3;  {extra length of unfolding intermediate in nm}
      EaUnf=28;   {unfolding activation energy in pN*nm}
var residual:real;
begin
dN:=Afreq*dt*N*exp(-EaUnf+F*dXunf)/kT;
i:=trunc(dN);
N:=N-i; {decrease number of native domains}
residual:=dN-i;   {but if the increase is fractional}
if residual>random then dec(N); {then decrement N on probablistic basis}
if N<0 then N:=0;
end;

procedure refold; {increases the number of folded domains}
const Afreq=1E8;  {attempts per second}
      EaFold=0;   {delta G for refolding in pN*nm}
      dXref=8;    {nm length decrease to make one fold}
var residual:real;
begin
dN:=Afreq*dt*(Nmax-N)*exp(-EaFold-F*dXref)/kT;
i:=trunc(dN);
N:=N+i;     {increase folded number}
residual:=dN-i;   {but if the increase is fractional}
if residual > random then inc(N); {then increment N on probablistic basis}
if N>Nmax then N:=Nmax;
end;

begin {main}
initializeScreen;
N:=Nmax;             {all domains folded to start}
L1:=N*natLength;     {contour length of native titin}
L2:=L-N*(unfLength); {contour length of denatured region}
pipetteX:=L1;        {starting pipette position}
trapX:=0;            {starting position of bead in trap}
direction:=1;        {begin stretch half-cycle}
F:=0;                {initial force = zero}
endOfCycle:=false;
ext:=pipetteX-trapX;
repeat
  repeat
  movePipette;
     repeat
     adjustForce;
     adjustTrap;
     until abs(trapAdjustment) < 1; {nanometers}
  unfold;
  refold;
  L1:=N*natLength;  L2:=L-N*(unfLength); {adjust contour lengths}
  i:=round(xMax*ext/L); j:=yMax-round(yMax*F/fMax);
  if N > oldN then begin setColor(green);circle(i,j,2); end;
  if N < oldN then begin setColor(red); circle(i,j,2); end;
  if N = oldN then putPixel(i,j,blue);
  oldN:=N;
  until endOfCycle;
  ch:=readkey;
  if ch='e' then initializeScreen; {erase screen}
  endOfCycle:=false;
until ch='q';                 {quit}
closeGraph;
end.