#!/bin/sh
#

# calculates macroscopic binding constants from a chelator titration using 
# monte-carlo type fitting procedure
# finds optimal fit by randomly choosing values of the parameters within an
# n-dimensional cube. the cube shrinks until the error sum improves less than
# mindiff
 
if test $# = 0
        then
                echo "usage:   kjell inputfile logk1_est logk2_est tries denominator"
                exit
fi
 
# the inputfile can look as follows:
# kdchel 2.493e-5	chelator kd in M
# cchel 28.8    chelator concentration in uM
# cprot 25.3    protein concentration in uM
# naapo 0.375	start value for absorption in apo chelator
# naholo 0.048	start value for absorption in holo chelator
# c_0 1         initial ion concentration in uM
# v_0 2500      volume in uL
# 0 5           added vol in uL and concentration in mM (for each of the following rows)
# 0.5379        absorbance
# 5 2.35        added vol in uL and concentration in mM (for each of the following rows)
# 0.4514        absorbance
# 0.3678        absorbance
# 0.2915        absorbance
#   .
#   .
#   .
# logk1_est and logk2_est are starting values for the macroscopic
# binding constants in 10logs (in log(M-1)).
# tries is the number of random values used
# denominator is the shrinkage factor of the cube
# anders malmendal 1997-2001
                                                                                                         
nawk 'BEGIN{
  k1_0=10^'$2'
  k2_0=10^'$3'
  k12_0=k1_0*k2_0
  fp_0=1.0	# protein concentration correction
  tries='$4'
  denom='$5'
  mindiff=0.0001
  k1_range=0.3
  k12_range=0.1
  fp_range=0.0
  aapo_range=0.1
  aholo_range=0.1
  pi=3.1415926359
  srand()
  minsum=100000
  print "round\tlog k1\t\tlog k2\t\tcorrection\taapo\t\taholo\t\terrsum" 
}

$1~/kdchel/ {kdchel=$2}
$1~/cprot/ {cprot_0=$2/1e6}
$1~/cchel/ {cchel_0=$2/1e6}
$1~/naapo/ || $1~/namax/ {aapo_0=$2}
$1~/naholo/ || $1~/namin/ {aholo_0=$2}
$1~/v_0/ {v_0=$2/1e6}
$1~/c_0/ {c_0=$2/1e6}

$1~/^[0-9.]/ && NF==1 {
a[++n]=$1
ionvol[n]=ionvol_
ioncons[n]=ioncons_
}

$1~/^[0-9.]/ && NF==2 {
  ionvol_=$1/1e6
  ioncons_=$2/1e3
}

END {
  nmax=n
  ctot=c_0
  vtot_=v_0
  nion=c_0*v_0
  for (n=1;n<=nmax;n++) {
    vtot_=vtot_+ionvol[n]
    vtot[n]=vtot_
    if (v_0!=0) korr[n] = vtot[n]/v_0
    if (korr[n]!=0)cchel[n]=cchel_0/korr[n]
    nion=nion+ioncons[n]*ionvol[n]
    iontot[n]=nion/vtot_
  }
  aapo_0=a[1]
  aholo_0=a[nmax]*korr[nmax]
  while (sqrt((1-oldsum/minsum)^2)>mindiff || oldsum==0) {
    oldsum=minsum
    if (++round!=1) {
      k1_best=k1_0
      k12_best=k12_0
      fp_best=fp_0
      aapo_best=aapo_0
      aholo_best=aholo_0 
    }
    for (i=1;i<=tries;i++) {
      sum=0
      k1=k1_0^(1-k1_range/2+k1_range*rand())
      k12=k12_0^(1-k12_range/2+k12_range*rand())
      k2=k12/k1
      fp=fp_0*(1-fp_range/2+fp_range*rand())
      aapo=aapo_0*(1-aapo_range/2+aapo_range*rand())
      aholo=aholo_0*(1-aapo_0/aholo_0*aholo_range/2+aapo_0/aholo_0*aholo_range*rand())
      y = 1e-12
      f2=2*k2
      for (n=1;n<=nmax;n++) {
        cprot = cprot_0*fp/korr[n]
        f1=cprot*k1
        slut=0 
        for (k=1;k<=100 && slut==0;k++) {
          t1=1/(y+kdchel)
          t2=1+f2*y 
          t3=1+k1*y*(1+k2*y)
          t4=1/t3
          par=(1+cchel[n]*t1+f1*t2*t4) 
          fy=iontot[n]-y*par
          dfy=-par-y*(f1*t4^2*(f2*t3-t2*(k1+f2*k1*y))-cchel[n]*t1^2)
          yn=y-fy/dfy
          if(sqrt((1-yn/y)^2)<1e-4) slut=1 
          y=yn
        } 
        acalc=(aapo-(aapo-aholo)*y/(y+kdchel))/korr[n]
        sum=sum+(a[n]-acalc)^2
      }
      if (sum<minsum) {
        minsum=sum
        k1_0=sqrt(k1*k1_0)
        k12_0=sqrt(k12*k12_0)
        fp_0=(fp+fp_0)/2
        aapo_0=(aapo+aapo_0)/2
        aholo_0=(aholo+aholo_0)/2
        k1_best=k1
        k12_best=k12
        fp_best=fp
        aapo_best=aapo
        aholo_best=aholo
#      print i,log(k1)/log(10),log(k2)/log(10),fp,aapo,aholo,sum
      }
    }
    if ( round%1==0 ) {
      print round "\t" log(k1_best)/log(10) "\t" log(k12_best/k1_best)/log(10) "\t" fp_best "\t" aapo_best "\t" aholo_best "\t" minsum
    }
    k1_0=k1_best
    k12_0=k12_best
    fp_0=fp_best
    aapo_0=aapo_best
    aholo_0=aholo_best
    k1_range=k1_range/denom
    k12_range=k12_range/denom
    fp_range=fp_range/denom
    aapo_range=aapo_range/denom
    aholo_range=aholo_range/denom
  }
}' $1 > $6