Executable code

File name: thermo.c

#include <math.h>
 
#include "bases.h"
#include "tables.h"
 
#define R (0.001986) /* Boltzmann's constant in kcal/(mol*K) */
 
/* Given the forward and reverse first-order (or pseudo-first-order)
   rate constants for a chemical reaction that is changing a reactant
   into a product, compute the fraction of initial reactants that will
   still exist after a given time from the start of the reaction.
 
   In this program, the "reactant" is a single empty or filled site on
   the template DNA strand, and the "product" is that same site in the
   opposite state (filled or empty). */
 
double compute_fraction_remaining (kf, kr, Dt)
     double kf;                 /* Forward first-order rate constant. */
     double kr;                 /* Reverse first-order rate constant. */
     double Dt;                 /* Length of time interval. */
{
  double kt = kf + kr;                /* Total rate constant. */
  return (kr + kf*exp(-kt*Dt)) / kt;  /* Fraction remaining. */
}
 
 
/* Given the standard change in Gibbs free energy for a reaction and
   given the temperature in Kelvins, compute the equilibrium constant
   for the reaction. */
 
double compute_keq (DG0, T)
     double DG0, T;
{
  double k = exp(-DG0/(R*T));
  return k;
}
 
/* Given oligos a and b, length n, the first 5'->3' and the second
   3'->5', and the temperature, estimate the equilibrium constant of
   binding of the two when aligned as given (with the 5' end of "a"
   binding to the 3' end of "b", and vice-versa, and everything in
   between aligned accordingly). */
 
double compute_DG0 (a, b, n, T, show_detail) 
     int *a, *b;
     int n;
     double T;
     int show_detail;
{
  double DH0,DS0,DG0;
  int i;
 
  if (n==0) return 0;
 
  /* Compute parameters for the first pair.
     This assumes that there is a GC-pair in the oligo, that
     Breslauer's DGi is correct, and that the enthalpy of initial
     pair formation is 0. */
  DH0 = 0;
  DS0 = -.01677;
 
  if (show_detail) {
 
    printf("         Contribution of doublet or mismatch");
    printf("        Running totals\n");
    printf("         ---------------------------------- ");
    printf("   ---------------------------\n");
    printf("   base   enthalpy   entropy   free energy  ");
    printf("   enthalpy  entropy   free energy\n");
    printf("   pair   kcal/mol  kcal/mol*K  kcal/mol    ");
    printf("  kcal/mol  kcal/mol*K  kcal/mol\n");
 
    printf("  ----- ---------- ---------- ----------    ");
    printf("---------- ---------- ----------\n");
    printf("                                            ");
    printf("%10g %10g %10g\n\n", DH0,DS0,DH0-T*DS0);
  }
 
  for (i=0; i<n; i++) {
    char mchar;
    double pDH0=DH0;
    double pDS0=DS0;
    int doublet=0;
    if (i<n-1 && (a[i] ^ b[i]) == 2 && (a[i+1] ^ b[i+1]) == 2) {
      doublet=1;
      /* We're on the first base-pair of a matching doublet, so add in
         the doublet parameters. */
 
      DH0 += doublet_enthalpy[a[i]][a[i+1]];
      DS0 += doublet_entropy[a[i]][a[i+1]];
 
    }
    if ((a[i] ^ b[i]) != 2) {
      mchar='x';
      /* We're on a mismatched base, so add in the parameters for the
         mismatched virtual stack. */
 
      DH0 += mismatch_enthalpy[a[i]][b[i]];
      DS0 += mismatch_entropy[a[i]][b[i]];
 
      /* Correction for context not being AxT/TyA. */
 
      if (i > 0) {
 
        /* Subtract out guessed contribution if there had been a T:A
           pair on left. */
 
        DH0 -= (doublet_enthalpy[base_T][a[i]] +
                doublet_enthalpy[b[i]][base_A]) / 2.0;
        DS0 -= (doublet_entropy[base_T][a[i]] +
                doublet_entropy[b[i]][base_A]) / 2.0;
 
        /* Add in guessed contribution from the actual pair on the left. */
 
        DH0 += (doublet_enthalpy[a[i-1]][a[i]] +
                doublet_enthalpy[b[i]][b[i-1]]) / 2.0;
        DS0 += (doublet_entropy[a[i-1]][a[i]] +
                doublet_entropy[b[i]][b[i-1]]) / 2.0;
 
      }
 
      if (i < n-1) {
 
        /* Subtract contribution from A:T on right. */
 
        DH0 -= (doublet_enthalpy[a[i]][base_A] +
                doublet_enthalpy[base_T][b[i]]) / 2.0;
        DS0 -= (doublet_entropy[a[i]][base_A] +
                doublet_entropy[base_T][b[i]]) / 2.0;
 
        /* Add contrib. for actual pair on right. */
 
        DH0 -= (doublet_enthalpy[a[i]][a[i+1]] +
                doublet_enthalpy[b[i+1]][b[i]]) / 2.0;
        DS0 -= (doublet_entropy[a[i]][a[i+1]] +
                doublet_entropy[b[i+1]][b[i]]) / 2.0;
 
      }
    } else mchar = '-';
    if (show_detail) {
      printf("%3d %c%c%c ",i+1,base_names[a[i]],mchar,base_names[b[i]]);
      if (doublet == 1) {
        printf("\n");
        printf("        %10g %10g %10g    %10g %10g %10g\n",DH0-pDH0,DS0-pDS0,
               (DH0-pDH0)-T*(DS0-pDS0),DH0,DS0,DH0-T*DS0);
      } else if (mchar == 'x') {
        printf("%10g %10g %10g    %10g %10g %10g\n",DH0-pDH0,DS0-pDS0,
               (DH0-pDH0)-T*(DS0-pDS0),DH0,DS0,DH0-T*DS0);
        printf("\n");
      } else {
        printf("\n\n");
      }
    }
  }
  /* Note: the above is an oversimplification of how mismatches at ends
     work. */
 
  DG0 = DH0 - T*DS0;
  return DG0;
}