particles in a box Master Equation

1 collaborator

Luis Mayorga (Author)

WHAT IS IT?

Simulation of a Reversible Reaction

This model simulates molecules in a thermal bath, where they switch between two states based on energy constraints. The states differ in standard enthalpy (H°) and entropy (S°), with:
- Height difference representing ΔH° (enthalpy change).
- Sector width representing ΔS° (entropy change).

This two-sector box provides a visual representation of a single-molecule reaction. In this simulation, a master equation is implemented to analytically estimate the proportion of molecules in each state as the reaction progresses from different initial conditions.

The master equation output enables the estimation of free energy across time, tracking its evolution toward equilibrium. This approach provides insights into the relationship between Probabilities (used within the Master Equation) and reaction dynamics, equilibrium behavior, and energy minimization principles.

HOW IT WORKS

Master Equation: an analytical prediction from State Transitions Probabilities

A molecule can be in State 1 (Box Left) or State 2 (Box Right), switching randomly between them.
The master equation calculate the number of molecules at time t+1 from the number of molecules in each state at time t, and the probabilities of switching from one state to the other. Knowing the initial number of molecules in each state, the master equation can be solved to predict the progresion of the reaction.

In the simulation, transitions probabilities depend on whether molecules acquire enough energy to overcome Ea, along with the probability that they will land in one compartment or the other.
Thus, state transitions are influenced by:
- Enthalpy difference (ΔH°) between states.
- Entropy difference (ΔS°), affecting molecular distribution.
- Energy of Activation (Ea), acting as a transition barrier. The notation should be Ebarr. The actual Ea is Ebarr or Ebarr + ΔH° (see below)

For a molecule in Box Left:
- If enthalpy1 ≤ enthalpy2, molecules must overcome Ea + ΔH° to move right.
- If enthalpy1 > enthalpy2, molecules require Ea to move right, but Ea - ΔH° to return left. Notice that in this case - ΔH° is a positive number

These probabilities are calculated using:

ifelse enthalpy1 <= enthalpy2  
  [ set alfa entropy2 * exp(- (Ea + deltaH) / R / temperature)  
    set beta entropy1 * exp(- (Ea) / R / temperature) ]  
  [ set alfa entropy2 * exp(- (Ea) / R / temperature)  
    set beta entropy1 * exp(- (Ea - deltaH) / R / temperature) ]

Here:
- alfa represents the probability of transitioning to State 2 (Box Right).
- beta represents the probability of transitioning to State 1 (Box Left).

MASTER EQUATION SOLUTION

Knowing the transition probabilities alfa and beta, the master equation can be solved, and the proportion of molecules in each state calculated as the reaction progresses to reach equilibrioum.

set P1 1 / (alfa + beta) * (beta + rr * exp(-(alfa + beta) * time ))  
set P2 1 / (alfa + beta) * (alfa - rr * exp(-(alfa + beta) * time ))

Where:
- P1 = Fraction of molecules in Box Left (State 1).
- P2 = Fraction of molecules in Box Right (State 2).
- rr = Adjustment for initial reaction conditions (reaction advance). - time (in tick number).

At equilibrium, P1 and P2 stabilize.

ESTIMATING FREE ENERGY FROM THE MASTER EQUATION

The system's free energy is dynamically estimated:

set averageEnergy R * temperature * ln (P2 / (P1 + 1E-10)) + deltaG0

Here:
- Boltzmann factor (ln(P2/P1)) quantifies state distribution.
- R * T accounts for thermal energy effects.
- ΔG0 represents standard free energy at equilibrium.

This ensures the model automatically tracks ΔG changes, showing free energy evolution toward equilibrium (ΔG → 0).

THINGS TO NOTICE

The master equation estimates equilibrium analytically, while stochastic simulation fluctuates around these predictions. Notice how fluctuations decrease when more molecules are present.

THINGS TO TRY

Modify Ea and observe delays
Adjust temperature and track free energy evolution
Alter enthalpy and entropy values
Change the extent of reaction

SCREENING PARAMETERS AND REPEATING SIMULATIONS RUNS

In the “Tools” tab of NetLogo, you will find “BehaviorSpace,” which enables you to explore different values of model parameters and perform repeated simulation runs.

KEY TAKEAWAYS

THE POWER OF THE MASTER EQUATION IN SYSTEM ANALYSIS

The master equation is a powerful tool for studying system behavior when transition probabilities are known. It allows a precise mathematical representation of dynamic changes, offering valuable insights into equilibrium, reaction kinetics, and free energy evolution.

COPYRIGHT AND LICENSE

CC BY-NC-SA 3.0

This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 License. To view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/3.0/ or send a letter to Creative Commons, 559 Nathan Abbott Way, Stanford, California 94305, USA.

Commercial licenses are also available. To inquire about commercial licenses, please contact Uri Wilensky at uri@northwestern.edu.

This model was created as part of the project: CONNECTED MATHEMATICS: MAKING SENSE OF COMPLEX PHENOMENA THROUGH BUILDING OBJECT-BASED PARALLEL MODELS (OBPML). The project gratefully acknowledges the support of the National Science Foundation (Applications of Advanced Technologies Program) -- grant numbers RED #9552950 and REC #9632612.

This model was converted to NetLogo as part of the projects: PARTICIPATORY SIMULATIONS: NETWORK-BASED DESIGN FOR SYSTEMS LEARNING IN CLASSROOMS and/or INTEGRATED SIMULATION AND MODELING ENVIRONMENT. The project gratefully acknowledges the support of the National Science Foundation (REPP & ROLE programs) -- grant numbers REC #9814682 and REC-0126227. Converted from StarLogoT to NetLogo, 2002.

Comments and Questions

Please start the discussion about this model! (You'll first need to log in.)

Click to Run Model

globals
[
  tick-delta                          ;; how much we advance the tick counter this time through
  max-tick-delta                      ;; the largest tick-delta is allowed to be
  avg-speed-init avg-energy-init      ;; initial averages
  avg-speed avg-energy                ;; current averages
  fast medium slow lost-particles     ;; current counts
  percent-medium percent-slow         ;; percentage of current counts
  percent-fast percent-lost-particles ;; percentage of current counts
  box-edge1
  box-edge2
  leftbox-edge
  rightbox
  floor1
  floor2
  entropy1
  ;entropy2
 ; temperature
  prev-temp
  R
  yboltz
  enthalpy1
 ; enthalpy2
 ; scale
  Eafloor
 ; initialstate ; particles in initially in state 1
  deltaG
  deltaG0
  deltaG0e
  -deltaG0e
  deltaS
  deltaH
  Kin1
  Kin2
  K
  tracer
  a
  initialState
  P1 ; probabilidad estado 1
  P2 ; probabilidad estado 2
  alfa ; probabilidad de saltar a estado 2 estando en 1
  beta ; probabilidad de saltar a estado 1 estando en 2
  rr ; factor probabilidades
  pinicial; probabilidad inicial de estar en estado 1
  time
  averageEnergy
  guessGreen
  guessRed
  uTotal

]

breed [ particles particle ]
breed [ flashes flash ]

flashes-own [birthday]

particles-own
[
  speed mass energy          ;; particle info
  last-collision
  state
]

to setup
  RESET-TICKS
  clear-all
  set-default-shape particles "circle"
  set-default-shape flashes "plane"
  set R 1.9858775 / 1000 ;constante de gases en kcal/mol/K
  set initialState 1 - avance
  set prev-temp temperature
 ; set entropy1 0.3
  set entropy1 1 - entropy2
  set enthalpy1 0
  ifelse enthalpy1 <= enthalpy2
   [set floor1 0
     set floor2 enthalpy2 * scale] ; scale 1kcal= scale y coordenate
   [set floor2 0
     set floor1 -1 * scale * enthalpy2] ; scale 1kcal=scale y coordinate
  set Eafloor max (list floor1 floor2) + Ea * scale
 ; set initialstate 0.1
  set box-edge2 (max-pxcor)
  set box-edge1 min-pxcor + (entropy1 * max-pxcor * 2)
  set deltaH enthalpy2 - enthalpy1
  set deltaS  R * ln (entropy2 / entropy1) ;calculo entropía
  set deltaG0 deltaH - temperature * deltaS ;calculo delta G cero teórica
  ; establece cantiadades para la ecuación maestra
  ifelse enthalpy1 <= enthalpy2
  ; si deltaH es positivo
  [set alfa entropy2 * exp(- (Ea + deltaH) / R / temperature);si deltaH es positivo, para ir a rojo tiene que saltar
  set beta entropy1 * exp(- (Ea) / R / temperature)
  ]
  ;si deltaH es negativo
  [set alfa entropy2 * exp(- (Ea) / R / temperature)
  set beta entropy1 * exp(- (Ea - deltaH) / R / temperature); si deltaH es negativo para ir a verde tiene que saltar
  ]
  set pinicial initialState
  set rr alfa * pinicial - beta * (1 - pinicial)
  make-box
  make-particles
  update-variables
  reset-ticks
  ;set time ticks

  set a 1
  set time 1
end 

to go
  ask particles [ bounce ]
    update-variables
    update-plots
    display
    set a a + 1
  set time time + 1
  if single = true [stop]
;  if a mod 300 = 0 [stop]
end 

to update-variables
  set deltaG0 deltaH - temperature * deltaS ;calculo delta G cero teórica
    if temperature != prev-temp [
    masterEquation
  set prev-temp temperature
  set time 0
  ]
  ; ecuación maestra
  set P1 1 / (alfa + beta) * (beta + rr * exp(-(alfa + beta) * time ))
  set P2 1 / (alfa + beta) * (alfa - rr * exp(-(alfa + beta) * time ))
  set averageEnergy R * temperature * ln (P2 / (P1 + 1E-10)) + deltaG0
  set Kin1 count turtles with [state = 1]
  set Kin2 count turtles with [state = 2]
  set K Kin2 / (Kin1 + 0.0000000000000001)

    ; potenciales químicos
  let uGreen  -1 * deltaG0 + R * temperature * ln (Kin1 + 0.0000001)
  let uRed R * temperature * ln (Kin2 + 0.01)
  set uTotal uGreen * Kin1 + uRed * Kin2

  set deltaG0e (- R * temperature * ln (K + 1E-10))
  set -deltaG0e -1 * deltaG0e
  set deltaG R * temperature * ln (K + 1E-10) + deltaG0
  ; ecuación maestra
  ;set guessGreen Kin1 + Kin2 * beta - Kin1 * alfa
  ;set guessRed Kin2 + Kin1 * alfa - Kin2 * beta
end 

to masterEquation
   ; establece cantiadades para la ecuación maestra
  ifelse enthalpy1 <= enthalpy2
  ; si deltaH es positivo
  [set alfa entropy2 * exp(- (Ea + deltaH) / R / temperature);si deltaH es positivo, para ir a rojo tiene que saltar
  set beta entropy1 * exp(- (Ea) / R / temperature)
  ]
  ;si deltaH es negativo
  [set alfa entropy2 * exp(- (Ea) / R / temperature)
  set beta entropy1 * exp(- (Ea - deltaH) / R / temperature); si deltaH es negativo para ir a verde tiene que saltar
  ]
  set Kin1 count turtles with [state = 1]
  set Kin2 count turtles with [state = 2]
  set pinicial Kin1 / (Kin1 + Kin2)
  set rr alfa * pinicial - beta * (1 - pinicial)
end 

to bounce  ;; particle procedure

 ifelse state = 1
     [
     set yboltz floor1 - scale * ln (random-float 1) * R * temperature
     if yboltz >= Eafloor
     [set xcor min-pxcor + random (max-pxcor - min-pxcor)]
     ]
     [
     set yboltz floor2 - scale * ln (random-float 1) * R * temperature
     if yboltz >= Eafloor
     [set xcor min-pxcor + random (max-pxcor - min-pxcor)]
     ]

  ifelse yboltz > max-pycor
       [set ycor  max-pycor]
       [set ycor yboltz]

  ifelse (xcor >= min-pxcor and xcor < box-edge1)
   [set state 1
    set color green]
  [set state 2
    set color red]
end 

;;;
;;; drawing procedures

to make-box

     ; white the part of the  box that is inactive
     ifelse floor1 <= floor2
     [
       ask patches with [ (pxcor > box-edge1) and (pycor < floor2) ]
     [ set pcolor white ]
     ]
     [
       ask patches with [ (pxcor < box-edge1) and (pycor < floor1) ]
     [ set pcolor white ]
     ]

       ; limite superior
    ask patches with [ pycor > max-pycor - 5 ]
    [ set pcolor blue ]

    ; limite inferior
    ask patches with [ (pycor < floor1 + 5 and pycor > floor1 ) and (pxcor <= box-edge1) ]
    [ set pcolor yellow ]

    ask patches with [ (pycor < floor2 + 5 and pycor > floor2) and (pxcor >= box-edge1) ]
    [ set pcolor yellow ]
    ; limite izquierdo
    ask patches with [ (pxcor < min-pxcor + 5) and (pycor >= floor1) ]
    [ set pcolor yellow ]

    ; limite derecho
    ask patches with [ (pxcor > max-pxcor - 5) and (pycor >= floor2) ]
    [ set pcolor yellow ]

    ; limite central
    ask patches with [ (pxcor < box-edge1 + 3 and pxcor > box-edge1 - 3) and (pycor <= Eafloor) ]
    [ set pcolor yellow ]
end 
;;;


;; creates initial particles

to make-particles
  create-particles number-of-particles
  [

    set size 6
    ifelse random-float 1 < initialstate
    [
      set xcor min-pxcor + random (box-edge1 - min-pxcor)
      set yboltz floor1 - scale * ln (random-float 1) * R * temperature
      ifelse yboltz > max-pycor
       [set ycor  max-pycor]
       [set ycor yboltz]
    ]

    [
      set xcor box-edge1 + random (box-edge2 - box-edge1)
      set yboltz floor2 - scale * ln (random-float 1) * R * temperature
      ifelse yboltz > max-pycor
       [set ycor  max-pycor]
       [set ycor yboltz]
    ]


   ifelse (xcor >= min-pxcor and xcor < box-edge1)
  [set state 1
    set color green]
  [set state 2
    set color red]
    ;;random-position

  ]
   set K count turtles with [state = 2] / (count turtles with [state = 1] + 0.00000000000000000000001)
  set deltaG0e (- R * temperature * ln (K + 0.0000000000001))
  set deltaG R * temperature * ln (K + 0.0000000000001) + deltaG0
end 




; Copyright 1997 Uri Wilensky.
; See Info tab for full copyright and license.

There are 11 versions of this model.

Uploaded by	When	Description	Download
Luis Mayorga	about 4 hours ago	new reaction	Download this version
Luis Mayorga	2 months ago	prueba	Download this version
Luis Mayorga	2 months ago	prueba	Download this version
Luis Mayorga	2 months ago	prueba	Download this version
Luis Mayorga	2 months ago	prueba1	Download this version
Luis Mayorga	2 months ago	Reverted to older version	Download this version
Luis Mayorga	2 months ago	Reverted to older version	Download this version
Luis Mayorga	2 months ago	To fit with a future paper	Download this version
Luis Mayorga	3 months ago	To fit with a future paper	Download this version
Luis Mayorga	3 months ago	prueba	Download this version
Luis Mayorga	3 months ago	Initial upload	Download this version

Attached files

File	Type	Description	Last updated
MasterEquation - Copy.nlogo	extension	MasterEq does not work in web (???) but it can be run in the computer	about 4 hours ago, by Luis Mayorga	Download
MasterEquation.nlogo	extension	MasterEq does not work in web	3 months ago, by Luis Mayorga	Download
particles in a box Master Equation.png	preview	Preview for 'particles in a box Master Equation'	3 months ago, by Luis Mayorga	Download

This model does not have any ancestors.

This model does not have any descendants.

NetLogo