/* compile with R CMD SHLIB protocols_propensities.c or (on Windows) Rcmd SHLIB protocols_propensities.c */

#include <R.h>
#include <Rmath.h>

#define Gy 2.

typedef enum {
  G_PTEN, Gi_PTEN, PTENt, PTEN,
  PIP, PIPp,
  Akt, Aktp,
  G_MDM2, Gi_MDM2, MDM2t, MDM2, MDM2p, MDM2pn,
  G_P53, Gi_P53, P53t, P53n, P53pn,
  IKKKn, IKKKa,
  IKKn, IKKa, IKKi, IKKii,
  G_A20, Gi_A20, A20t, A20,
  G_IkBa, Gi_IkBa, IkBat, IkBa, IkBap, IkBan, 
  NFkB, NFkBn,
  NFkB_IkBa, NFkB_IkBap, NFkB_IkBan,
  TNF, Gamma,
  Rec, Reci,
  DSB,
  AF
} enum_places;


/************************************
*             Protocol
*/
double TNF_before_IR(double dStartTime, double *y) {
  dStartTime = fround(dStartTime,0);
  if (dStartTime == 0) {
    y[Gamma] = 0;
    y[TNF] = 0;
    return 30*3600;
  } else if (dStartTime == 30*3600) {
    y[TNF] = 10;
    return 33*3600;
  } else if (dStartTime == 33*3600) {
    y[Gamma] = Gy;
    return 34*3600;
  } else if (dStartTime == 34*3600) {
    y[Gamma] = 0;
    y[TNF] = 0;
  }

  return 82*3600;
}


/************************************
*             Propensities
*/
#define d0 3e-5
#define d1 1.5e-4
#define h0 7.

double MDM2_degr(double time, double *y) {
  double DSB_sq = y[DSB]*y[DSB];
  return (d0 + d1 * DSB_sq / (h0*h0 + DSB_sq )) * y[MDM2];
}

double MDM2p_degr(double time, double *y) {
  double DSB_sq = y[DSB]*y[DSB];
  return (d0 + d1 * DSB_sq / (h0*h0 + DSB_sq )) * y[MDM2p];
}

double MDM2pn_degr(double time, double *y) {
  double DSB_sq = y[DSB]*y[DSB];
  return (d0 + d1 * DSB_sq / (h0*h0 + DSB_sq )) * y[MDM2pn];
}

#define a0 1e-4
#define a1 1e-3

double P53n_phos(double time, double *y) {
  double DSB_sq = y[DSB]*y[DSB];
  return (a0 + a1 * DSB_sq / (h0*h0 + DSB_sq )) * y[P53n];
}

#define d3 1e-4
#define d4 1e-13

double P53n_degr(double time, double *y) {
  return (d3 + d4*y[MDM2pn]*y[MDM2pn]) * y[P53n];
}

#define d5 1e-4
#define d6 1e-14

double P53pn_degr(double time, double *y) {
  return (d5 + d6*y[MDM2pn]*y[MDM2pn]) * y[P53pn];
}

#define ka 1e-4
#define ka20 1e4

double IKKKn_activ(double time, double *y) {
  //  return ka * ka20 / (ka20 + y[A20]) * y[Rec] * y[IKKKn];
  // As ka * ka20 == 1...
  return y[Rec] * y[IKKKn] / (ka20 + y[A20]);
}

#define k2 1e4
#define k3 3e-3

double IKKa_inactiv(double time, double *y) {
  return k3 * y[IKKa] * (k2 + y[A20]) / k2;
}

#define nc1 1e-1
#define h101 6e4

double A20t_transcr(double time, double *y) {
  return nc1 * h101 / (h101 + y[P53pn]) * y[G_A20];
}

double IkBat_transcr(double time, double *y) {
  return nc1 * h101 / (h101 + y[P53pn]) * y[G_IkBa];
}

#define p1 200
#define q3 2.5e-11
#define q4 1

double AF_synth(double time, double *y) {
  double tmp = q3 * y[P53pn]*y[P53pn];
  return p1 * tmp/(q4 + tmp);
}

#define q0 1e-4
#define q1 5e-13

double Gi_MDM2_activ(double time, double *y) {
  return (q0 + q1 * y[P53pn]*y[P53pn]) * y[Gi_MDM2];
}

double Gi_PTEN_activ(double time, double *y) {
  return (q0 + q1 * y[P53pn]*y[P53pn]) * y[Gi_PTEN];
}

#define a6 1e-2

#define Th 0.65

#define d9_over_p1 1e-6

double DSB_increase_AF(double time, double *y) {
  double sgn;
  double tmp = y[AF] * d9_over_p1 - Th;
  if (tmp > 0)
    sgn = 1;
  else if (tmp < 0)
    sgn = -1;
  else
    sgn = 0;

  return a6 * (sgn + 1);
}

#define Nsat 3.
#define dRep 2e-14

double DSB_decrease(double time, double *y) {
  return dRep * y[P53pn]*y[P53pn] * y[DSB]/(y[DSB] + Nsat);
}