#lang racket
(require (planet williams/science/ode-initval)
         plot racket/gui plot/utils)

; Simulation of the enzyme reaction>
;   S + E <-> ES -> P
; S  substrate
; E  enzyme
; ES enzyme-substrate-complex (below the letter c is used)
; P  product

; The model
;   s' = - k1 * s * e
;   e' = - k1 * s * e + k_minus1 * c + k2 * c 
;   c' =   k1 * s * e - k_minus1 * c - k2 * c
;   p' =   k2 * c
; Here k1      is the reaction constant of S+E -> ES
;      kminus1 is the reaction constant of S+E <- ES
;      k2      is the reaction constant of        ES -> P

; The simulation begins with the following concentrations:
(define smax 100.0) ; TODO find better way
(define emax 10.0)
(define cmax 10.0)
(define pmax 10.0)

(define steps 1000)
(define time-end 1.5)
(define h 0.001)

(define s0 30.0)
(define e0 1.0)
(define c0 0.0)
(define p0 0.0)
; The parameters are:
(define k1 10.0)
(define kminus1 100)
(define k2 1.0)
; And the theoretical values for vmax and km are:
(define vmax (* e0 k2))
(define km (/ (+ kminus1 k2) k1))

(define (func t y f params)
  (let ((k1      (first params))
        (kminus1 (second params))
        (k2      (third params))
        (s (vector-ref y 0))
        (e (vector-ref y 1))
        (c (vector-ref y 2))
        (p (vector-ref y 3)))
    (vector-set! f 0 (->inexact (- (* k1 s e))))
    (vector-set! f 1 (->inexact (+ (- (* k1 s e))    (* kminus1 c)     (* k2 c))))
    (vector-set! f 2 (->inexact (+    (* k1 s e)  (- (* kminus1 c)) (- (* k2 c)))))
    (vector-set! f 3 (->inexact                                       (* k2 c)))))

(define (->inexact x)
  (if (inexact? x)
      x 
      (exact->inexact x)))

(define frame  (new frame% [label "Michaelis-Menten"]))
(new text-field% 
     [parent frame] 
     [label "k1"]   
     [init-value (number->string k1)]
     [callback (λ (t e) (let ([x (string->number (send t get-value))])
                          (when x (set! k1 x) (draw))))])
(new text-field% 
     [parent frame] [label "k_-1"] [init-value (number->string kminus1)]
     [callback (λ (t e) (let ([x (string->number (send t get-value))])
                          (when x (set! kminus1 x) (draw))))])
(new text-field% [parent frame] [label "k2"] [init-value (number->string kminus1)]
     [callback (λ (t e) (let ([x (string->number (send t get-value))])
                          (when x (set! k2 x) (draw))))])
(new text-field% [parent frame] [label "tidslut"] [init-value (number->string time-end)]
     [callback (λ (t e) (let ([x (string->number (send t get-value))])
                          (when x (set! time-end x) (draw))))])

(define pane   (new vertical-pane% [parent frame] [stretchable-height #f]))
(define panel2 (new horizontal-pane% [parent frame]))
(define canvas  (new canvas% [parent panel2] [min-width 400] [min-height 400] [paint-callback (λ (c dc) (draw))]))
(define canvas2 (new canvas% [parent panel2] [min-width 400] [min-height 400] [paint-callback (λ (c dc) (draw))]))

(define update-sliders values) ; make-slider will update this value

(define (make-slider parent the-label min max the-callback get-value get-max-value)
  (let ([pane (new horizontal-pane% [parent parent])])
    (new message% [parent pane] [label the-label] [min-width 20])
    (new slider% 
         [parent pane]
         [label ""]
         [style '(horizontal plain)]
         [init-value (inexact->exact (floor (* max (/ (get-value) (get-max-value)))))]
         [min-value min]
         [max-value max]
         [min-width 100]
         [callback the-callback])
    (let ([m (new message% 
                  [parent pane] 
                  [label (number->string (get-value))]
                  [min-width 20])])
      (set! update-sliders
            (compose (λ (x) 
                       (send m set-label (number->string (get-value))))
                     update-sliders)))))

(define slider-max 10000)

(define substrate-slider 
  (make-slider pane
               "[S]" 0 slider-max 
               (lambda (slider event)
                 (when (eq? event 'slider)
                   (set! s0 (exact->inexact (* smax (/ slider-max) (send slider get-value)))))
                 (draw))
               (λ () s0)
               (λ () smax)))

(define enzyme-slider 
  (make-slider pane
               "[E]" 0 slider-max 
               (lambda (slider event)
                 (when (eq? event 'slider)
                   (set! e0 (exact->inexact (* emax (/ slider-max) (send slider get-value)))))
                 (draw))
               (λ () e0)
               (λ () emax)))

(define complex-slider 
  (make-slider pane
               "[C]" 0 slider-max 
               (lambda (slider event)
                 (when (eq? event 'slider)
                   (set! c0 (exact->inexact (* cmax (/ slider-max) (send slider get-value)))))
                 (draw))
               (λ () c0)
               (λ () cmax)))

(define product-slider 
  (make-slider pane
               "[P]" 0 slider-max 
               (lambda (slider event)
                 (when (eq? event 'slider)
                   (set! p0 (exact->inexact (* pmax (/ slider-max) (send slider get-value)))))
                 (draw))
               (λ () p0)
               (λ () pmax)))

(define (draw)
  (let* ((type rk4-ode-type)
         (dim 4)
         (step (make-ode-step type dim))
         (system (make-ode-system func #f dim (list k1 kminus1 k2)))
         (t 0.0)
         (y (vector s0 e0 c0 p0))
         (y-err (make-vector dim))
         (dydt-in (make-vector dim))
         (dydt-out (make-vector dim))
         (s-values '())
         (e-values '())
         (c-values '())
         (p-values '()))
    (ode-system-function-eval system t y dydt-in)
    (let loop ()
      (when (< t time-end)
        (ode-step-apply step t h
                        y y-err
                        dydt-in
                        dydt-out
                        system)
        (unless (or (ormap nan? (vector->list y))
                    (ormap infinite? (vector->list y)))
          (set! s-values (cons (vector t (vector-ref y 0)) s-values))
          (set! e-values (cons (vector t (vector-ref y 1)) e-values))
          (set! c-values (cons (vector t (vector-ref y 2)) c-values))
          (set! p-values (cons (vector t (vector-ref y 3)) p-values)))
        (vector-set! dydt-in 0 (vector-ref dydt-out 0))
        (vector-set! dydt-in 1 (vector-ref dydt-out 1))
        (vector-set! dydt-in 2 (vector-ref dydt-out 2))        
        (vector-set! dydt-in 3 (vector-ref dydt-out 3))
        (set! t (+ t h))
        (loop)))
    (update-sliders 'update!)
    (plot/dc (list (points s-values #:color 0 #:size 3 #:label "[S]")
                               (points e-values #:color 1 #:size 2 #:label "[E]")
                               (points c-values #:color 2 #:size 2 #:label "[ES]")
                               (points p-values #:color 3 #:size 2 #:label "[P]"))
             (send canvas get-dc)
             0 0
             (send canvas get-width)
             (send canvas get-height)
             #:x-min 0.0
             #:x-max time-end
             #:x-label "Tid"
             #:y-min 0.0
             #:y-max (* 1.5 (apply max (map (λ (v) (vector-ref v 1)) (append e-values c-values p-values))))
             #:y-label "Koncentration")
    (let* ([ss (map (λ (v) (vector-ref v 1)) s-values)]
           [ps (map (λ (v) (vector-ref v 1)) p-values)]
           [vs (map (λ (x) (/ x h)) (map - (cons (car ps) ps) (append ps (list (car (reverse ps))))))]
           [pvs (map vector ss (rest vs))])
      (plot/dc (list (points pvs #:color 0 #:size 2 #:label "v")
                     (function (λ (s) (/ (* vmax s) (+ km s)))))
               (let ([dc (send canvas2 get-dc)])
                 (send dc set-font (make-object font% 50 'swiss)) 
                 dc)
               0 0
               (send canvas get-width)
               (send canvas get-height)
               #:x-min (apply min ss)
               #:x-max (* 1.5 (apply max ss))
               #:x-label "Substratkoncentration  [S]"
               #:y-min (apply min vs)
               #:y-max (if (= (apply min vs) (apply max vs))
                           (+ (apply min vs) 1.0)
                           (* 1.5 (apply max vs)))
               #:y-label "Hastighed  v(t)"
               #:title "Produktkoncentrationens vækstrate som funktion af substratkoncentrationen"))))

(send frame show #t)
(draw)