#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)