Team:Technion Israel/Modeling
3-\(\alpha \)-HSD Kinetic Model
\(\begin{array}{l}\frac{{d\left[ {3 - \alpha - d} \right]}}{{dt}} = \frac{{d{{\left[ {3 - \alpha - d} \right]}_{forward}}}}{{dt}} + \frac{{d{{\left[ {3 - \alpha - d} \right]}_{backward}}}}{{dt}} = \\ = \left[ {3 - \alpha - HSD} \right] \cdot \left( {{K_{ca{t_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_{{m_1}}} + \left[ {DHT} \right]}} - {K_{ca{t_2}}} \cdot \frac{{\left[ {3 - \alpha - d} \right]}}{{{K_{{m_2}}} + \left[ {3 - \alpha - d} \right]}}} \right)\end{array}\)
Background
3-\(\alpha \)-HSD is a name of a group of enzymes which convert certain hormones (like DHT) to another hormone (like 3-α-diol) and vice versa by means of oxidation and reduction. There are several strands of this enzyme, with different level of potency. In humans, the enzyme is encoded by the AKR1C4 gene, while in rat it is encoded by the AKR1C9 gene. We chose the rat version of the enzyme because it is more efficient in breaking down DHT [need article].
All AKRs catalyze an ordered bi-bi reaction in which the cofactor binds first, followed by the binding of the steroid substrate. The steroid product is the first to leave, and the cofactor is the last. In this mechanism, \({K_{cat}}\) represents the slowest step in the kinetic sequence [2].
Approaches to modeling the process
1. Cofactor saturation assumption
We assume that the levels of the cofactors on the scalp are high enough so that they are always saturated in the enzymatic reaction. The advantage of this approach is that we can use the michaelis-menten reversible equation to describe the reaction. As we will explain later, this assumption may not be correct so we will offer other approaches. Another major disadvantage of this model is that it does not take the levels of cofactors into consideration, so it cannot help us to predict the system's behavior for different cofactor concentrations.
2. New Model Development
Taking cofactors into consideration, we can use principles from statistical mechanics in order to develop a completely new enzyme reaction function. The advantage of this approach is that it describes the kinetics of the enzyme in much more detail than michaelis-menten reversible, and can even offer some explanations to our wetlab results. The disadvantage of this approach is that there are no reaction constants available for it, so we will have to estimate them.
Approach 1 – cofactor saturation assumption
If we assume that the levels of cofactor in the enzyme's environment are high enough that they become saturated, the probability of finding an enzyme that is not connected to a cofactor is negligible. We also need to assume that the concentrations of both cofactors are almost equal, so the inhibition effect [need article] will not affect the reaction (as we will show later, a large ratio in favor of one cofactor will inhibit the other direction of the reaction). The new kinetic schematic is:
Since both levels of NADPH and NADP are saturated, we'll assume product inhibition occurs only with DHT and \(3\alpha diol\), so the reaction will resemble michaelis-menten reversible. Since the degradation rates of the hormones on the scalp are unknown, we will neglect them by assuming the degradation is slower by several orders of magnitude than the enzymatic reaction.
We can summarize the reactions by the following coupled differential equations:
\[\left( 1 \right)\left\{ {\begin{array}{*{20}{c}}{\frac{{d\left[ {3\alpha diol} \right]}}{{dt}} = \frac{{{V_{{m_f}}} \cdot \frac{{\left[ {DHT} \right]}}{{{k_s}}} - {V_{{m_r}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{k_p}}}}}{{1 + \frac{{\left[ {DHT} \right]}}{{{k_s}}} + \frac{{\left[ {3\alpha diol} \right]}}{{{k_p}}}}}}\\{\frac{{d\left[ {DHT} \right]}}{{dt}} = \frac{{{V_{{m_r}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{k_p}}} - {V_{{m_f}}} \cdot \frac{{\left[ {DHT} \right]}}{{{k_s}}}}}{{1 + \frac{{\left[ {DHT} \right]}}{{{k_s}}} + \frac{{\left[ {3\alpha diol} \right]}}{{{k_p}}}}}}\end{array}} \right.\]
Where:
- \({V_{{m_f}}}\) is the maximum forward reaction rate. Attained when all enzyme molecules are bound to the substrate (DHT).
- \({V_{{m_r}}}\) is the maximum backward reaction rate. Attained when all enzyme molecules are bound to the product(\(3\alpha Diol\)).
- \({K_s}\) is the substrate concentration at which the forward reaction rate is at half-maximum.
- \({K_p}\) is the product concentration at which the backward reaction rate is at half-maximum.
While we couldn't find the kinetic constants for human scalp, we found an article which measured them on rat skin [5]. We will assume the constants are on the same order of magnitude as on the human scalp.
Time domain simulation
We simulated the system described above. The simulation has been done using the following constants:
| Parameter | Value | Units | source | comment |
|---|---|---|---|---|
| \({V_{{m_f}}}\) | 5.63 | \(\frac{{nmol}}{{\min }}\) | Calculated from [5] | For \({10_{mg}}\) of enzyme |
| \({V_{{m_r}}}\) | 16.28 | \(\frac{{nmol}}{{\min }}\) | Calculated from [5] | For \({10_{mg}}\) of enzyme |
| \({K_s}\) | 0.38 | \(\mu M\) | From [5] | |
| \({K_p}\) | 2.79 | \(\mu M\) | From [5] |
In order to understand the breakdown process of DHT by the enzyme, we simulated the system for three initial concentrations of the hormone:
Graphs from top to bottom are the simulation for the following initial concentrations of DHT:
\(0.1 \cdot {K_s}\), \(1 \cdot {K_s}\) and \(10 \cdot {K_s}\)
From the simulation we can see that for an initial concentration that is low by an order of magnitude from \({K_s}\), the time for the system to reach stable state is almost the same as for an initial concentration of \({K_s}\). However for a value that is higher by an order of magnitude, the time increases significantly.
The following video shows a parameter scan for different initial concentration values:
From this video we can learn that the time to break down initial concentrations that are larger by an order of magnitude from \({K_s}\) increase significantly for every small increment. By looking at equation (1) we can find two possible reasons for this:
- DHT saturation - for every increment in the initial concentration, the value of the derivative approaches \({V_{{{\mathop{\rm m}\nolimits} _f}}}\). For values that are larger by an order of magnitude than \({K_s}\), large increments in the initial concentration do not change much the value of the derivative.
- Product inhibition - most of the DHT molecules are converted to \(3\alpha - diol\). Larger concentration of the product decreases the breakdown rate of the substrate.
We can use this information in order to determine the necessary enzyme concentration required to break down an initial concentration of DHT at a certain time (Increasing the enzyme level will increase \({V_{\max }}\)).
Percentage of DHT breakdown
Let C be the initial concentration of DHT, and the initial concentration of \(3\alpha - diol\). The relation between the substrate and the product is:
(2) \(\left[ {DHT} \right] + \left[ {3\alpha - diol} \right] = C\)
At stable state, \(\frac{d}{{dt}} = 0\), so from equation (1) we get: \[\begin{array}{l}\frac{{{V_{{m_f}}} \cdot \frac{{{{\left[ {DHT} \right]}_{final}}}}{{{k_s}}} - {V_{{m_r}}} \cdot \frac{{{{\left[ {3\alpha diol} \right]}_{final}}}}{{{k_p}}}}}{{1 + \frac{{{{\left[ {DHT} \right]}_{final}}}}{{{k_s}}} + \frac{{{{\left[ {3\alpha diol} \right]}_{final}}}}{{{k_p}}}}} = 0\\{V_{{m_f}}} \cdot \frac{{{{\left[ {DHT} \right]}_{final}}}}{{{k_s}}} - {V_{{m_r}}} \cdot \frac{{{{\left[ {3\alpha diol} \right]}_{final}}}}{{{k_p}}}\mathop = \limits^{(2)} {V_{{m_f}}} \cdot \frac{{{{\left[ {DHT} \right]}_{final}}}}{{{k_s}}} - {V_{{m_r}}} \cdot \frac{{C - {{\left[ {DHT} \right]}_{final}}}}{{{k_p}}} = 0\\{\left[ {DHT} \right]_{final}} \cdot \left( {\frac{{{V_{{m_f}}}}}{{{k_s}}} + \frac{{{V_{{m_r}}}}}{{{k_p}}}} \right) = \frac{{{V_{{m_r}}} \cdot C}}{{{k_p}}}\\{\left[ {DHT} \right]_{final}} = C \cdot \frac{{{V_{{m_r}}}}}{{\frac{{{k_p}}}{{{k_s}}} \cdot {V_{{m_f}}} + {V_{{m_r}}}}}\\ \Rightarrow \left[ \% \right] = \left( {1 - \frac{{{{\left[ {DHT} \right]}_{final}}}}{C}} \right) \cdot 100 = \left( {1 - \frac{{{V_{{m_r}}}}}{{\frac{{{k_p}}}{{{k_s}}} \cdot {V_{{m_f}}} + {V_{{m_r}}}}}} \right) \cdot 100\mathop = \limits^{from\,\,\,table} 0.717\% \end{array}\]
Problems with the approach
Looking at figure 1, we can see that there is an element of cofactor inhibition of the enzyme that this model does not take into account. In order for the enzyme to convert DHT to \(3\alpha - diol\) (or vice versa), the specific cofactor of the reaction has to be connected to the enzyme, transforming it so it can bind only to the substrate or product. As we will see later, in order for a michaelis menten reversible reaction to occur, both cofactor levels need to be high and at a certain ratio. The fact that this assumption may not be correct is one of the main reasons we chose to Overproduce NADPH as a part of our project.
In fact, this model does not treat the fact that the enzyme requires a cofactor in order to work, so we cannot use it as a good simulator of our system.
Nevertheless, in situations where these conditions exist, this model offers a simple way for prediction that is easy to understand and analyse.
Approach 2 – cofactor saturation assumption
In this part we will look thoroughly at the mechanics of the enzyme and develop new rate equations for the system that will take cofactor levels into account. Later we will analyse those equations and explain how they correlate with our previous model and with our wetlab results.
The rate equation for an enzyme
Let \({K_{cat}}\) be the the number of substrate molecule each enzyme site converts to product per unit time and P the number of product molecules. If we have N enzymes in our solution, the rate equation for the product of the enzyme will be:
\(\frac{{dP}}{{dt}} = {K_{cat}} \cdot \sum\limits_{i = 1}^N {{1_{enzym{e_i}}}} \) where:
\({1_{enzym{e_i}}} = \left\{ {\begin{array}{*{20}{c}}{1,}&{if\,the\,enzme\,currently\,converts\,a\,substrate}\\{0,}&{else}\end{array}} \right.\)
is the indicator fuction for each enzyme.
If an enzyme operates in both directions, with different velocities for each direction, the rate equation will be:
\[\frac{{dP}}{{dt}} = \sum\limits_{i = 1}^N {{K_{enzym{e_i}}}} \] where:
\[{K_{enzym{e_i}}} = \left\{ {\begin{array}{*{20}{c}}{{K_{ca{t_{forward}}}}}&{if\,enzyme\,converts\,substrate\,to\,product}\\{{K_{ca{t_{backwards}}}}}&{if\,enzyme\,converts\,product\,to\,substrate}\\0&{else}\end{array}} \right.\]
We'll assume that \({P_{reactio{n_{forward}}}}\) and \({P_{reactio{n_{backwards}}}}\) - the probabilities for each reaction to occur - are the same for all the enzymes in the system. As we will see later, They are a function of a number of factors, among them is the level of the substrate in the solution and the probability of binding a substrate molecule to the enzyme.
Since there is usually a very large number of molecules in the solution, we can use the law of large numbers[put a random article about it here]:
\[\begin{array}{l}\,\,\,\,\,\,\,\,\,\frac{{dP}}{{dt}} = \sum\limits_{i = 1}^N {{K_{enzym{e_i}}}} \approx N \cdot E\left[ {{K_{enzym{e_i}}}} \right] = N \cdot \left( {{K_{ca{t_{forward}}}} \cdot {P_{reactio{n_{forward}}}} - {K_{ca{t_{backward}}}} \cdot {P_{reactio{n_{backward}}}}} \right)\\ \Rightarrow \frac{{d\left[ P \right]}}{{dt}} = \left[ {enzyme} \right] \cdot \left( {{K_{ca{t_{forward}}}} \cdot {P_{reactio{n_{forward}}}} - {K_{ca{t_{backward}}}} \cdot {P_{reactio{n_{backward}}}}} \right) = {V_{{m_f}}} \cdot {P_{reactio{n_{forward}}}} - {V_{{m_r}}} \cdot {P_{reactio{n_{backward}}}}\end{array}\]
Next, we will find \({P_{reactio{n_{forward}}}}\) and \({P_{reactio{n_{backwards}}}}\) for our enzyme.
Thermodynamic model
The AKR1C9 has one site for binding a cofactor. Once the cofactor is bound, the enzyme has a site which the appropriate hormone binds to. Let us look at a system with one receptor that has two binding sites and four ligands:
We'll assume we have \({L_A}\) ligands from type A, \({L_B}\) ligands from type B and so on for C and D. The solution has \(\Omega \) "volume cells" which can contain only one ligand.
To find the probability for the enzyme to be in a certain state, we will use the Boltzmann distribution[insert book explanation pages]\[({e^{ - \frac{{{\varepsilon _{macrostate}}}}{{KT}}}})\], which is a probability distribution that gives the probability that a system will be in a certain state as a function of that state’s energy and the temperature of the system.
The probability for the system to be in a certain macrostate is:
\[{P_{macrostat{e_i}}} = \frac{{\overbrace {{e^{ - \frac{{{\varepsilon _{macrostat{e_i}}}}}{{KT}}}}}^{{\rm{The}}\,{\rm{Weight}}\,{\rm{of}}\,{\rm{the}}\,{\rm{macrostate}}}}}{{\underbrace {\sum\limits_{j = 1}^M {{e^{ - \frac{{{\varepsilon _{macrostat{e_j}}}}}{{KT}}}}} }_{{\rm{The}}\,{\rm{sum}}\,{\rm{of}}\,{\rm{all}}\,{\rm{macrostates}}\,{\rm{weights}}}}}\]Where M is the number of states accessible to the system.
The weight for each macrostate is the sum of the weights of all of its microstates:
\[{e^{ - \frac{{{\varepsilon _{macrostate}}}}{{KT}}}} = \sum\limits_i {{e^{ - \frac{{{\varepsilon _{microstat{e_i}}}}}{{KT}}}}} \]A microstate of the system is one possible arrangement of \(\left( {{L_A} + {L_B} + {L_C} + {L_D}} \right)\) ligands in \[\Omega \] cells or/and the binding sites of the receptor. Every macrostate of the system has a certain number of these arrangements.
Each ligand can have two stable states of energy:
- When the ligand is free in the solution - \({\varepsilon _{sol}}\). For the sake of simplicity we will assume all four ligands have the same energy when free in the solution.
- When the ligand is bound to a site in the receptor - \({\varepsilon _b}\). Each ligand will have different energy from the others when bound to the enzyme.
The sum of the energies of all the ligands will give us the energy of a microstate:
\({\varepsilon _{microstate}} = \sum\limits_{i = 1}^{{L_A}} {{\varepsilon _{{A_i}}}} + \sum\limits_{i = 1}^{{L_B}} {{\varepsilon _{{B_i}}}} + \sum\limits_{i = 1}^{{L_C}} {{\varepsilon _{{C_i}}}} + \sum\limits_{i = 1}^{{L_D}} {{\varepsilon _{{D_i}}}} \)In our system, only two ligands that are bound to the enzyme can have the energy \({\varepsilon _b}\). We assume the rest have the energy \({\varepsilon _{sol}}\), no matter in which box they are in. For that reason, we can say that all the microstates of each macrostate are the same.
For example, for the macrostate in which only one ligand of type A is bound to the enzyme:
\(\begin{array}{l}{\varepsilon _{microstate}} = {\varepsilon _{sol}} \cdot \left[ {\left( {{L_A} - 1} \right) + {L_B} + {L_C} + {L_D}} \right] + {\varepsilon _{{b_A}}}\\ \Rightarrow {e^{ - \frac{{{\varepsilon _{macrostate}}}}{{KT}}}} = \sum\limits_{i = 1}^{MP} {{e^{ - \frac{{{\varepsilon _{microstat{e_i}}}}}{{KT}}}}} = \sum\limits_{i = 1}^{MP} {{e^{ - \frac{{{\varepsilon _{microstate}}}}{{KT}}}}} = MP \cdot {e^{ - \frac{{{\varepsilon _{microstate}}}}{{KT}}}} = MP \cdot {e^{ - \frac{{{\varepsilon _{sol}} \cdot \left[ {\left( {{L_A} - 1} \right) + {L_B} + {L_C} + {L_D}} \right] + {\varepsilon _{{b_A}}}}}{{KT}}}}\end{array}\)Where MP is the multiplicity of the microstate – the number of possible arrangements of ligands in the system.
Counting microstates
We'll notice that for every arrangement of type A ligands, there are numerous possible arrangements of the other ligands. Also, there is redundancy in the total possible arrangements – since all ligands of type A are the same, it does not matter which one occupies which cell. Same thing for particles of types B, C and D.
For example, the number of microstates of the system for the macrostate in which none of the ligands are bound to the receptor is:
\[\begin{array}{l}\left( {\begin{array}{*{20}{c}}{\Omega - {L_B} - {L_C} - {L_D}}\\{{L_A}}\end{array}} \right) \cdot \left( {\begin{array}{*{20}{c}}{\Omega - {L_A} - {L_C} - {L_D}}\\{{L_B}}\end{array}} \right) \cdot \left( {\begin{array}{*{20}{c}}{\Omega - {L_A} - {L_B} - {L_D}}\\{{L_C}}\end{array}} \right) \cdot \left( {\begin{array}{*{20}{c}}{\Omega - {L_A} - {L_B} - {L_C}}\\{{L_D}}\end{array}} \right) = \\ = \frac{{\left( {\Omega - {L_B} - {L_C} - {L_D}} \right)!}}{{{L_A}!\, \cdot \,\left( {\Omega - {L_B} - {L_C} - {L_D} - {L_A}} \right)!}} \cdot \frac{{\left( {\Omega - {L_A} - {L_C} - {L_D}} \right)!}}{{{L_B}!\, \cdot \,\left( {\Omega - {L_A} - {L_C} - {L_D} - {L_B}} \right)!}} \cdot \frac{{\left( {\Omega - {L_A} - {L_B} - {L_D}} \right)!}}{{{L_C}!\, \cdot \,\left( {\Omega - {L_A} - {L_B} - {L_D} - {L_C}} \right)!}} \cdot \frac{{\left( {\Omega - {L_A} - {L_B} - {L_C}} \right)!}}{{{L_D}!\, \cdot \,\left( {\Omega - {L_A} - {L_B} - {L_C} - {L_D}} \right)!}}\mathop \approx \limits^{\left( 1 \right)} \\ \approx \frac{{{{\left( {\Omega - {L_B} - {L_C} - {L_D}} \right)}^{{L_A}}} \cdot {{\left( {\Omega - {L_A} - {L_C} - {L_D}} \right)}^{{L_B}}} \cdot {{\left( {\Omega - {L_A} - {L_B} - {L_D}} \right)}^{{L_C}}} \cdot {{\left( {\Omega - {L_A} - {L_B} - {L_C}} \right)}^{{L_D}}}}}{{{L_A}!\, \cdot \,{L_B}! \cdot {L_C}! \cdot {L_D}!}}\end{array}\]Let us look at all the possible macrostates of the system: