Executable Code

File name: sim.c

#include <stdio.h>
#include <math.h>               /* for sqrt() */
#include "struct.h" /* for oligo struct */
#include "thermo.h" /* for compute_k */
#include "bases.h" /* for print_bases (temporary) */
#include "print.h"  /* print_parse */
#include "sim.h"                /* REAL_TIME, MEDIUM_SPEED, FASTEST */
 
/* frand() returns doubles in [0,1) */
#define frand() ((double)random()/(32768.0*65536.0))
 
static double tprof(T,t)
     double T;                  /* Goal (annealing) temperature. */
     double t;                  /* Time since start of cycle. */
{
  double aT;
  if (t < 60.0) {
    aT = (371.15 * (1.0 - (t/60.0))) + (T * (t/60.0));
  } else {
    aT = T;
  }
  return aT;
}
 
#define QUIET 0
#define VERBOSE 1
 
static void bprof (rstrand, len, oligos, nOligos, syms, nsyms, T, f,
                   verbose)
     int *rstrand;
     int len;
     Oligo *oligos;
     int nOligos;
     Sym *syms;
     int nsyms;
     double T,f;
     int verbose;
{
  int *done;
  double top_p,best_p,k,p,best_k,best_DG0,best_lkf;
  int posi,oli,best_posi,best_oli,i;
  int hit_low=0;
 
  printf("\n");
  printf("   ");
  print_parse(rstrand,len,REVERSE_ORDER,syms,nsyms);
  printf("\n");
  printf("3'-");
  print_bases(rstrand,len);
  printf("-5'\n");
 
  done = (int *)malloc(len*nOligos*sizeof(int));
  top_p = 0.0;
  for(posi=0;posi<len;posi++){
    for(oli=0;oli<nOligos;oli++){
      done[posi+oli*len] = 0;
    }
  }
  for(;;){
    best_p = 0.0;
    for(posi=0;posi<len;posi++){
      for(oli=0;oli<nOligos;oli++){
        if (!done[posi+oli*len]) {
          if (posi+oligos[oli].len <= len) {
            double lkf;
            double DG0 = compute_DG0(oligos[oli].bases, &rstrand[posi],
                                     oligos[oli].len, T, DONT_PRINT);
            k = compute_keq(DG0,T);
            p = k*f/(1+k*f);
            lkf = log10(k*f);
            if (p>best_p) {
              best_p = p; best_posi = posi; best_oli = oli;
              best_k = k; best_DG0 = DG0; best_lkf = lkf;
              if (best_p > top_p) top_p = best_p;
            }
          }
        }
      }
    }
    if (best_p < .01) {
      if (hit_low>=3) break;
      printf("\nBelow .01 probability of binding:\n");
      hit_low++;
    }
 
    printf("\n");
 
    /* Print matches and mismatches. */
    printf("   ");
    for(i=0;i<best_posi;i++) putchar(' ');
    for(i=0;i<oligos[best_oli].len;i++){
      if (rstrand[i+best_posi] != (oligos[best_oli].bases[i]^2))
        putchar('X');
      else
        putchar(':');
    }
    printf("    DG0 = %g\n",best_DG0);
 
    /* Print base sequence. */
    for(i=0;i<best_posi;i++) putchar(' ');
    printf("5'-");
    print_bases ( oligos[best_oli].bases, oligos[best_oli].len );
    printf("-3' k = %g\n",best_k);
 
    /* Print parse. */
    printf("   ");
    for(i=0;i<best_posi;i++) putchar(' ');
    print_parse (oligos[best_oli].bases, oligos[best_oli].len,
                 NORMAL_ORDER, syms, nsyms);
    printf("    p = %f\n",best_p);
 
    done[best_posi+best_oli*len] = 1;
 
    if (verbose) {
      /* Show the detailed base-by-base computation of DG0. */
      printf("\n");
      compute_DG0(oligos[best_oli].bases,&rstrand[best_posi],
                  oligos[best_oli].len, T, SHOW_DETAIL);
    }
 
  }
  free(done);
  printf("\n");
}
 
 
 
/* Take as input an initial strand, an array of oligos,
   and an annealing temperature, and an oligo concentration,
   and simulate away! */
 
void sim (istrand, len, oligos, nOligos, syms, nsyms, goalT, f, nticks,
          ncycles, speed)
     int *istrand;
     int len;
     Oligo *oligos;
     int nOligos;
     Sym *syms;
     int nsyms;
     double goalT;              /* Annealing temp, in Kelvin */
     double f;                  /* Concentration of each oligo, in Molar */
     int speed;                 /* REAL_TIME, MEDIUM_SPEED, or FASTEST */
{
  int *strand;                  /* 5'->3' current strand */
  int *rstrand;                 /* 3'->5' reverse of current strand */
  int order = 0;                /* Print the template in 5'->3' order.
                                   1 would mean to print in 3'->5'. */
  int i;
  int tick;
  int *occupy;
  int cycle;
 
  strand = (int *)malloc(len*sizeof(int));
  rstrand = (int *)malloc(len*sizeof(int));
  occupy = (int *)malloc(len*sizeof(int));
 
  /* Copy the initial strand to the current strand. */
  for(i=0;i<len;i++) strand[i] = istrand[i];
 
  for(cycle=1;((ncycles==0) || (cycle<=ncycles));cycle++){
    int j,o; 
 
    /* Compute 3'->5' reverse of current strand. */
    for(i=0;i<len;i++)
      rstrand[i] = strand[(len-1)-i];
 
    if (ncycles==0) {
      char linebuf[100];
      linebuf[0] = 'i';
      while(linebuf[0]!='\n'){
        printf("\n");
        printf("Command: ");
        fgets(linebuf,99,stdin);
        switch(linebuf[0]){
        case 'h':
          printf("\n");
          printf("Commands:\n");
          printf("\n");
          printf("h        Help.\n");
          printf("?        Oligo binding profile.\n");
          printf("v        Same but with more detail.\n");
          printf("n        No animation.\n");
          printf("r        Real time.\n");
          printf("RETURN   Do one cycle.\n");
          break;
        case '?':
          bprof(rstrand,len,oligos,nOligos,syms,nsyms,goalT,f,QUIET);
          break;
        case 'v':
          bprof(rstrand,len,oligos,nOligos,syms,nsyms,goalT,f,VERBOSE);
          break;
        case 'n':
          speed = NO_ANIMATION;
          printf("No animation.\n");
          break;
        case 'r':
          speed = REAL_TIME;
          printf("Animate in real time.\n");
          break;
        case 'm':
          speed = MEDIUM_SPEED;
          printf("Animate at medium speed (approx. 30x reality).\n");
          break;
        case 's':
          speed = FASTEST;
          printf("Animate at fastest speed.\n");
          break;
        }
      }
      putchar('\n');
    }
 
    /* Initialize the occupation vector.  It contains -2 in unoccupied
       spaces, -1 in occupied spaces, except that the first base of
       an occupied space contains the  of the oligo occupying that
       space. */
    for(i=0;i<len;i++) occupy[i] = -2;
 
    /* Annealing step: let oligos & DNA stew, for a long time. */
    for(tick=0;tick<nticks;tick++){
      int oi,i,j;
      int olen;
      int stick;
      int subtick;
      double T = tprof(goalT,(double)tick);
      double DG0;               /* Change std. Gibbs free energy. */
      double ka;                /* pseudo-1st-order rate constant for ass'n */
      double kd;                /* ditto for disassociation */
      double p;                 /* Probability of staying in same state. */
 
      /*Do "Poisson updating" of all sites. */
      for(subtick=0;subtick<len;subtick++){
 
        /* Pick a random site. */
        i = random()%len;
 
        /* If there is an oligo already at that site, compute
           whether to remove it. */
 
        if (occupy[i] >= 0) {
 
          /* Identify the oligo at this location and its length. */
 
          oi = occupy[i];
          olen = oligos[oi].len;
 
          /* Compute the change in standard Gibbs free energy
             for formation of this duplex. */
 
          DG0 = compute_DG0 (oligos[oi].bases, &rstrand[i], olen, T,
                             DONT_PRINT );
 
          /* Compute the first-order rate constant for annealing for
             oligos of this length. */
 
          ka = f * 3.5e5 / sqrt((double)olen);
 
          /* Compute the dissociation rate from the association
             rate and the first-order equilibrium constant. */
 
          kd = ka / (f * compute_keq (DG0, T));
 
          /* Compute the fraction of sites like this one that would
             still be filled after one second of the dissociation
             reaction for this duplex. */
 
          p = compute_fraction_remaining (kd, ka, 1.0);
 
          /* If we're not within the fraction of sites that stay
             filled, then dissociate. */
 
          if (frand() >= p ) {
            /*printf("Oligo %d dissociates from position %d.\n",oi,j);*/
            for(j=i;j<i+olen;j++) occupy[j] = -2;
          }
 
        } else if (occupy[i] == -2) {
 
          /* Otherwise, if the site is empty, then consider
             annealing with each oligo (once each, on average). */
 
          for (stick=0; stick<nOligos && occupy[i]==-2; stick++) {
 
            oi = random()%nOligos; /* Pick a random oligo. */
            olen = oligos[oi].len; /* Get its length. */
 
            /* If site is free, compute whether oligo will stick. */
 
            if (i+olen<=len) {
              int free = 1;
              for(j=i;j<i+olen;j++)
                if (occupy[j] != -2) { free=0; break; }
              if (free) {
 
                /* Compute the change in standard Gibbs free energy
                   for formation of this duplex. */
 
                DG0 = compute_DG0 (oligos[oi].bases, &rstrand[i], olen, T,
                                   DONT_PRINT );
 
                /* Compute the first-order rate constant for annealing for
                   oligos of this length. */
 
                ka = f * 3.5e5 / sqrt((double)olen);
 
                /* Compute the dissociation rate from the association
                   rate and the first-order equilibrium constant. */
 
                kd = ka / (f * compute_keq (DG0, T));
 
                /* Compute the fraction of sites like this one that
                   would still be unfilled by this oligo after one
                   second of the association reaction for this
                   duplex. */
 
                p = compute_fraction_remaining (ka, kd, 1.0);
 
                /* If we're not one of the still-unfilled fraction,
                   then associate. */
 
                if (frand() >= p) {
                  /*printf("Oligo %d sticks at position %d.\n",oi,j);*/
                  occupy[i] = oi;
                  for(j=i+1;j<i+olen;j++) occupy[j] = -1;
                }
              }
            }
          }
        }
      } /* End of sticking loop. */
 
      if (speed == REAL_TIME) {
        usleep(1000000);
      } else if (speed == MEDIUM_SPEED) {
        usleep(33333);
      }
      if (speed != NO_ANIMATION) printf("\033[H\033[2J");
      if (speed != NO_ANIMATION || tick == nticks-1) {
        printf("Cycle = %d, time = %d:%d%d, temp = %6.2f.\n",cycle,
               (tick+1)/60,((tick+1)%60)/10,(tick+1)%10,T-273.15);
        print_binders (rstrand, len, occupy, oligos, syms, nsyms);
      }
 
    } /* End of ticking loop. */
 
    /* Incorporate oligos into the new strand. */
    for(i=0;i<len;i++){
      if (occupy[i] >= 0) {
        int q;
        o = occupy[i]; j=0;
        for(j=0;j<oligos[o].len;j++){
          if (strand[i+j] != oligos[o].bases[j]) {
            strand[i+j] = oligos[o].bases[j];
          }
        }
        i+=oligos[o].len-1;
      } else {
        strand[i] = (rstrand[i]^2);
      }
    }
  } /* End of PCR loop. */
  printf("\n");
} /* End of sim function. */