Team:Technion Israel/Modeling

Team: Technion 2015


3\(\alpha \)-HSD Kinetic Model


3\(\alpha \)-HSD is the name of a group of enzymes which convert certain hormones (like DHT) to another hormone (like \(3\alpha-diol\)) and vice versa, by means of oxidation and reduction. There are several strands of this enzyme, with different levels of potency. In humans, the enzyme is encoded by the AKR1C4 gene, while in rats it is encoded by the AKR1C9 gene. We chose the rat version of the enzyme because it is more efficient in breaking down DHT [7].
Figure 1 – the kinetic scheme for 3\(\alpha \)-HSD

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].

Figure 2 – A schematic from an article [4] of the kinetic mechanisms for AKRs.

Approaches to modeling the process

1. Cofactor saturation assumption

We assume that the levels of the cofactors on the scalp are high enough that they are always at saturation 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 as well. Another major disadvantage of this model is that it does not take the levels of the cofactors into consideration, so it cannot help us 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 enzymatic 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 for our wet-lab 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 are at saturation, 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 inhibitory effect [need article] will not affect the reaction (as we will show later, a large ratio of one cofactor in relation to another will inhibit the other direction of the reaction). The new kinetic schematic is:

Figure 3 – kinetic schematic for high cofactor concentrations

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 a Michaelis-Menten reversible reaction. 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( I \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.\]


  • \({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 relevant for the human scalp, we found an article which measured them on rat skin [5]. We will assume that the constants are of 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:

Figure 4 – DHT and \(3\alpha - diol\) concentrations [\(\mu M\)] vs. time [hours].

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 lower by an order of magnitude from \({K_s}\), the time elapsed for the system to reach steady-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 elapsed to break down initial concentrations that are larger by an order of magnitude from \({K_s}\) increases significantly for every small increment. By looking at equation (I) 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\). A 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:

\(\left[ {DHT} \right] + \left[ {3\alpha - diol} \right] = C\)

At stable state, \(\frac{d}{{dt}} = 0\), so from equation (I) 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 acknowledge 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 tool 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}}}} \)


\({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 are usually a very large number of molecules in the solution, we can use the law of large numbers:
\[\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.

Developing the model

The drivation of this model is based on a similar model for binding of ligands to receptors found in the book[8].

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 of the enzyme being in a certain state, we will use the Boltzmann distribution[as explained in [8] pages 219-237]:

\[{e^{ - \frac{{{\varepsilon _{macrostate}}}}{{KT}}}}\]

The Boltzmann distribution 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:

  1. None of the Ligands are bound to the receptor:
    • All the ligands have an energy of \({\varepsilon _{sol}}\).
    • Energy of the macrostate: \(\left( {{L_A} + {L_B} + {L_C} + {L_D}} \right) \cdot {\varepsilon _{sol}}\)
    • Multiplicity of the state (the number of microstates that the macrostate contains) is as mentioned above:
      \[M{P_1} = \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}!}}\]
    • Weight of the macrostate (= Energy*Multiplicity): \({W_1} = M{P_1} \cdot {e^{ - \beta \cdot \left( {{L_A} + {L_B} + {L_C} + {L_D}} \right) \cdot {\varepsilon _{sol}}}}\) where \(\beta = \frac{1}{{KT}}\).
    • Because later we will want to normalize the weight of this macrostate to 1, we will multiply each macrostate weight by \[\frac{1}{{{W_1}}} = \frac{{{e^{ + \beta \cdot \left( {{L_A} + {L_B} + {L_C} + {L_D}} \right) \cdot {\varepsilon _{sol}}}}}}{{M{P_1}}}\], so that \({W_{{N_1}}} = 1\).
  2. One of the ligands of type A (cofactor) is bound to the receptor:
    • One ligand of type A is bound to its designated site in the receptor, and thus have an energy of \({\varepsilon _{{b_A}}}\).
    • \({L_A} - 1\) ligands of type A are unbound and thus have an energy of \({\varepsilon _{sol}}\).
    • All other ligands are unbound and thus have an energy of \({\varepsilon _{sol}}\).

    \[\begin{array}{l} \Rightarrow {\rm{ The}}\,{\rm{Energy}}\,{\rm{of}}\,{\rm{the}}\,{\rm{state }} = \left( {{L_A} - 1 + {L_B} + {L_C} + {L_D}} \right) \cdot {\varepsilon _{sol}} + {\varepsilon _{{b_A}}}\\\\{\rm{The}}\,{\rm{Multiplicity}}\,{\rm{of}}\,{\rm{the}}\,{\rm{state: }}\\{\rm{M}}{{\rm{P}}_2} = \frac{{{{\left( {\Omega - {L_B} - {L_C} - {L_D}} \right)}^{{L_A} - 1}} \cdot {{\left( {\Omega - {L_A} - 1 - {L_C} - {L_D}} \right)}^{{L_B}}} \cdot {{\left( {\Omega - {L_A} - 1 - {L_B} - {L_D}} \right)}^{{L_C}}} \cdot {{\left( {\Omega - {L_A} - 1 - {L_B} - {L_C}} \right)}^{{L_D}}}}}{{({L_A} - 1)!\, \cdot \,{L_B}! \cdot {L_C}! \cdot {L_D}!}}\\ \Rightarrow {\rm{ The}}\,{\rm{Weight}}\,{\rm{of}}\,{\rm{the}}\,{\rm{state: }}{{\rm{W}}_2} = {\rm{M}}{{\rm{P}}_2} \cdot {e^{ - \beta \cdot \left[ {\left( {{L_A} - 1 + {L_B} + {L_C} + {L_D}} \right) \cdot {\varepsilon _{sol}} + {\varepsilon _{{b_A}}}} \right]}}\end{array}\]

    \[\begin{array}{l}{\rm{Normalize}}\,{\rm{the}}\,{\rm{weight:}}\\{W_{{N_2}}} = \frac{{{{\left( {\Omega - {L_B} - {L_C} - {L_D}} \right)}^{{L_A} - 1}} \cdot {{\left( {\Omega - {L_A} - 1 - {L_C} - {L_D}} \right)}^{{L_B}}} \cdot {{\left( {\Omega - {L_A} - 1 - {L_B} - {L_D}} \right)}^{{L_C}}} \cdot {{\left( {\Omega - {L_A} - 1 - {L_B} - {L_C}} \right)}^{{L_D}}}}}{{({L_A} - 1)!\, \cdot \,{L_B}! \cdot {L_C}! \cdot {L_D}!}} \cdot {e^{ - \beta \cdot \left[ {\left( {{L_A} - 1 + {L_B} + {L_C} + {L_D}} \right) \cdot {\varepsilon _{sol}} + {\varepsilon _{{b_A}}}} \right]}} \cdot \\ \cdot \frac{{{L_A}!\, \cdot \,{L_B}! \cdot {L_C}! \cdot {L_D}!}}{{{{\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}}}}} \cdot {e^{ + \beta \cdot \left( {{L_A} + {L_B} + {L_C} + {L_D}} \right) \cdot {\varepsilon _{sol}}}} = \\ = \frac{{{L_A}}}{{\Omega - {L_B} - {L_C} - {L_D}}} \cdot {\left( {\frac{{\overbrace {\Omega - {L_A} - {L_C} - {L_D}}^{ > > 1} - 1}}{{\Omega - {L_A} - {L_C} - {L_D}}}} \right)^{{L_B}}} \cdot {\left( {\frac{{\overbrace {\Omega - {L_A} - {L_B} - {L_D}}^{ > > 1} - 1}}{{\Omega - {L_A} - {L_B} - {L_D}}}} \right)^{{L_C}}} \cdot {\left( {\frac{{\overbrace {\Omega - {L_A} - {L_B} - {L_C}}^{ > > 1} - 1}}{{\Omega - {L_A} - {L_B} - {L_C}}}} \right)^{{L_D}}} \cdot {e^{ - \beta \cdot \left( {{\varepsilon _{{b_A}}} - {\varepsilon _{sol}}} \right)}} \approx \\ \approx \frac{{{L_A}}}{{\Omega - {L_B} - {L_C} - {L_D}}} \cdot {\left( {\overbrace {\frac{{\Omega - {L_A} - {L_C} - {L_D}}}{{\Omega - {L_A} - {L_C} - {L_D}}}}^1} \right)^{{L_B}}} \cdot {\left( {\frac{{\overbrace {\Omega - {L_A} - {L_B} - {L_D}}^1}}{{\Omega - {L_A} - {L_B} - {L_D}}}} \right)^{{L_C}}} \cdot {\left( {\frac{{\overbrace {\Omega - {L_A} - {L_B} - {L_C}}^1}}{{\Omega - {L_A} - {L_B} - {L_C}}}} \right)^{{L_D}}} \cdot {e^{ - \beta \cdot \left( {{\varepsilon _{{b_A}}} - {\varepsilon _{sol}}} \right)}} = \\ = \frac{{{L_A}}}{{\Omega - {L_B} - {L_C} - {L_D}}} \cdot {e^{ - \beta \cdot \left( {{\varepsilon _{{b_A}}} - {\varepsilon _{sol}}} \right)}}\mathop \approx \limits^{\left( 2 \right)} \frac{{{L_A}}}{\Omega } \cdot {e^{ - \beta \cdot \left( {{\varepsilon _{{b_A}}} - {\varepsilon _{sol}}} \right)}}\mathop = \limits^{\left( 3 \right)} \frac{{\left[ A \right]}}{{\left[ {{c_0}} \right]}} \cdot {e^{ - \beta \cdot \left( {{\varepsilon _{{b_A}}} - {\varepsilon _{sol}}} \right)}}\mathop = \limits^{\left( 4 \right)} \frac{{\left[ A \right]}}{{{K_A}}}\end{array}\]

  3. Two ligands are bound to the receptor – one of type A and one of type B (Forward enzyme reaction):
    • \({L_A} - 1\) ligands of type A are unbound and thus have an energy of \({\varepsilon _{sol}}\).
    • \({L_B} - 1\) ligands of type B are unbound and thus have an energy of \({\varepsilon _{sol}}\).
    • One ligand of type A is bound to its designated site in the receptor, and thus have an energy of \({\varepsilon _{{b_A}}}\).
    • One ligand of type B is bound to it's designated site in the receptor, and thus have an energy of \({\varepsilon _{{b_B}}}\).
    • All other ligands are unbound and thus have an energy of \({\varepsilon _{sol}}\).
    \[\begin{array}{l} \Rightarrow {\rm{ The}}\,{\rm{Energy}}\,{\rm{of}}\,{\rm{the}}\,{\rm{state }} = \left( {{L_A} + {L_B} + {L_C} + {L_D} - 2} \right) \cdot {\varepsilon _{sol}} + {\varepsilon _{{b_A}}} + {\varepsilon _{{b_B}}}\\\\{\rm{The}}\,{\rm{Multiplicity}}\,{\rm{of}}\,{\rm{the}}\,{\rm{state: }}\\{\rm{M}}{{\rm{P}}_3} = \frac{{{{\left( {\Omega - {L_B} - 1 - {L_C} - {L_D}} \right)}^{{L_A} - 1}} \cdot {{\left( {\Omega - {L_A} - 1 - {L_C} - {L_D}} \right)}^{{L_B} - 1}} \cdot {{\left( {\Omega - {L_A} - 1 - {L_B} - 1 - {L_D}} \right)}^{{L_C}}} \cdot {{\left( {\Omega - {L_A} - 1 - {L_B} - 1 - {L_C}} \right)}^{{L_D}}}}}{{({L_A} - 1)!\, \cdot \,({L_B} - 1)! \cdot {L_C}! \cdot {L_D}!}}\\ \Rightarrow {\rm{ The}}\,{\rm{Weight}}\,{\rm{of}}\,{\rm{the}}\,{\rm{state: }}{{\rm{W}}_3} = {\rm{M}}{{\rm{P}}_3} \cdot {e^{ - \beta \cdot \left[ {\left( {{L_A} + {L_B} + {L_C} + {L_D} - 2} \right) \cdot {\varepsilon _{sol}} + {\varepsilon _{{b_A}}} + {\varepsilon _{{b_B}}}} \right]}}\\ \Rightarrow {\rm{ The}}\,{\rm{normalized}}\,{\rm{weight}}\,{\rm{(using}}\,{\rm{the}}\,{\rm{same}}\,{\rm{methods}}\,{\rm{as}}\,{\rm{for}}\,{\rm{the}}\,{\rm{previous}}\,{\rm{state): }}{{\rm{W}}_{{N_3}}} = \frac{{\left[ A \right]}}{{{K_A}}} \cdot \frac{{\left[ B \right]}}{{{K_B}}}\end{array}\]
  4. One of the ligands of type C (cofactor) is bound to the receptor:
    • Using the same methods as above, the normalized weight will be: \({W_{{N_4}}} = \frac{{\left[ C \right]}}{{{K_C}}}\)

  5. Two ligands are bound to the receptor – one of type C and one of type D (backwards enzyme reaction):
    • Using the same methods as above, the normalized weight will be: \({W_{{N_5}}} = \frac{{\left[ C \right]}}{{{K_C}}} \cdot \frac{{\left[ D \right]}}{{{K_D}}}\)

Approximations and definitions used in the development process:

\[\begin{array}{l}\left( 1 \right)\frac{{\Omega !}}{{\left( {\Omega - L} \right)!}}\mathop \approx \limits^{\Omega \gg L} {\Omega ^L}\\\left( 2 \right)\Omega - {L_{A,B,C,D}}\mathop \approx \limits^{\Omega \gg L} \Omega \Rightarrow \frac{{{L_{A,B,C,D}}}}{{\Omega - {L_{A,B,C,D}}}} \approx \frac{{{L_{A,B,C,D}}}}{\Omega }\\\left( 3 \right)\Omega = \left[ {{c_0}} \right] \cdot V,{L_i} = \left[ i \right] \cdot V\\\left( 4 \right){K_i} \equiv \left[ {{c_0}} \right] \cdot {e^{ + \beta \cdot \left( {{\varepsilon _{{b_i}}} - {\varepsilon _{sol}}} \right)}}\end{array}\]

To sum it all up, we can use the weight functions we found to find \({P_{reactio{n_{forward}}}}\) and \({P_{reactio{n_{backwards}}}}\):

\[\begin{array}{l}P_{forward} = {P_3} = \frac{{{W_{{N_3}}}}}{{\sum\limits_{i = 1}^5 {{W_{{N_i}}}} }} = \frac{{\frac{{\left[ A \right]}}{{{K_A}}} \cdot \frac{{\left[ B \right]}}{{{K_B}}}}}{{1 + \frac{{\left[ A \right]}}{{{K_A}}} + \frac{{\left[ A \right]}}{{{K_A}}} \cdot \frac{{\left[ B \right]}}{{{K_B}}} + \frac{{\left[ C \right]}}{{{K_C}}} + \frac{{\left[ C \right]}}{{{K_C}}} \cdot \frac{{\left[ D \right]}}{{{K_D}}}}}\\{P_{reverse}} = {P_5} = \frac{{{W_{{N_5}}}}}{{\sum\limits_{i = 1}^5 {{W_{{N_i}}}} }} = \frac{{\frac{{\left[ C \right]}}{{{K_C}}} \cdot \frac{{\left[ D \right]}}{{{K_D}}}}}{{1 + \frac{{\left[ A \right]}}{{{K_A}}} + \frac{{\left[ A \right]}}{{{K_A}}} \cdot \frac{{\left[ B \right]}}{{{K_B}}} + \frac{{\left[ C \right]}}{{{K_C}}} + \frac{{\left[ C \right]}}{{{K_C}}} \cdot \frac{{\left[ D \right]}}{{{K_D}}}}}\end{array}\]

And the overall reaction rate:

\(\frac{{d\left[ P \right]}}{{dt}} = \frac{{{V_{{m_f}}} \cdot \frac{{\left[ A \right]}}{{{K_A}}} \cdot \frac{{\left[ B \right]}}{{{K_B}}} - {V_{{m_r}}} \cdot \frac{{\left[ C \right]}}{{{K_C}}} \cdot \frac{{\left[ D \right]}}{{{K_D}}}}}{{1 + \frac{{\left[ A \right]}}{{{K_A}}} + \frac{{\left[ A \right]}}{{{K_A}}} \cdot \frac{{\left[ B \right]}}{{{K_B}}} + \frac{{\left[ C \right]}}{{{K_C}}} + \frac{{\left[ C \right]}}{{{K_C}}} \cdot \frac{{\left[ D \right]}}{{{K_D}}}}}\)

Since this is our model for the enzyme, we will write the equation again with our parameters:

\(\frac{{d\left[ {3\alpha diol} \right]}}{{dt}} = \frac{{d\left[ {NADP} \right]}}{{dt}} = - \frac{{d\left[ {DHT} \right]}}{{dt}} = - \frac{{d\left[ {NADPH} \right]}}{{dt}} = \frac{{{V_{{m_f}}} \cdot \frac{{\left[ {NADPH} \right]}}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} - {V_{{m_r}}} \cdot \frac{{\left[ {NADP} \right]}}{{{K_2}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}}}}{{1 + \frac{{\left[ {NADPH} \right]}}{{{K_1}}} + \frac{{\left[ {NADPH} \right]}}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} + \frac{{\left[ {NADP} \right]}}{{{K_2}}} + \frac{{\left[ {NADP} \right]}}{{{K_2}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}}}}\)

Compatibility with known theories

Here are some of the reasons why we think our model is correct:

  1. Compatibility with single-direction Michaelis-Menten:

    For \(\left[ {NADP} \right] < < {K_2}\) and \(\left[ {3\alpha diol} \right] < < {K_P}\), we get the following rate equation:

    \[\begin{array}{l}\frac{{d\left[ P \right]}}{{dt}} \approx \frac{{{V_{{m_f}}} \cdot \frac{{\left[ {NADPH} \right]}}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}}}}{{1 + \frac{{\left[ {NADPH} \right]}}{{{K_1}}} + \frac{{\left[ {NADPH} \right]}}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}}}} = \\ = {V_{{m_f}}} \cdot \frac{{\frac{{\left[ {DHT} \right]}}{{{K_s}}}}}{{\frac{{{K_1}}}{{\left[ {NADPH} \right]}} + 1 + \frac{{\left[ {DHT} \right]}}{{{K_s}}}}}\mathop \approx \limits^{\left[ {NADPH} \right] > > {K_1}} {V_{{m_f}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s} + DHT}}\end{array}\]

    This means that for very low levels of \(3\alpha-diol\) and NADP, the reaction is similar to Michaelis-Menten. In reality, the reaction will increase the levels of both substrates, so we cannot assume the entire reaction is similar to Michaelis-Menten, rather only the beginning of it. As we'll see later, this correlates with our wet-lab results for varying concentrations of DHT.

  2. Compatibility with Michaelis-Menten reversible reaction:

    If we assume both cofactor levels are saturated, and are high enough so that the reaction of the enzyme does not change their level much, we can treat them as constants:

    \[\begin{array}{l}\frac{{{V_{{m_f}}} \cdot \frac{{\left[ {NADPH} \right]}}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} - {V_{{m_r}}} \cdot \frac{{\left[ {NADP} \right]}}{{{K_2}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}}}}{{1 + \frac{{\left[ {NADPH} \right]}}{{{K_1}}} + \frac{{\left[ {NADPH} \right]}}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} + \frac{{\left[ {NADP} \right]}}{{{K_2}}} + \frac{{\left[ {NADP} \right]}}{{{K_2}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}}}}\mathop \approx \limits^{\left[ {NADPH} \right] \approx [NADP] \approx const = C > > {K_{1,}}{K_2}} \\ \approx \frac{{{V_{{m_f}}} \cdot \frac{C}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} - {V_{{m_r}}} \cdot \frac{C}{{{K_2}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}}}}{{1 + \frac{C}{{{K_1}}} + \frac{C}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} + \frac{C}{{{K_2}}} + \frac{C}{{{K_2}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}}}} = \frac{{{V_{{m_f}}} \cdot \frac{1}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} - {V_{{m_r}}} \cdot \frac{1}{{{K_2}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}}}}{{\frac{1}{{{K_1}}} + \frac{1}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} + \frac{1}{{{K_2}}} + \frac{1}{{{K_2}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}}}} = \\ = \frac{{{V_{{m_f}}} \cdot \frac{1}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} - {V_{{m_r}}} \cdot \frac{1}{{{K_2}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}}}}{{\left( {\frac{1}{{{K_1}}} + \frac{1}{{{K_2}}}} \right) + \frac{1}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} + \frac{1}{{{K_2}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}}}} = \frac{{{V_{{m_f}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} \cdot \frac{{{K_2}}}{{{K_1} + {K_2}}} - {V_{{m_r}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}} \cdot \frac{{{K_1}}}{{{K_1} + {K_2}}}}}{{1 + \frac{{\left[ {DHT} \right]}}{{{K_s}}} \cdot \frac{{{K_2}}}{{{K_1} + {K_2}}} + \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}} \cdot \frac{{{K_1}}}{{{K_1} + {K_2}}}}} = \\ = \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}'}}}}\end{array}\]

    This is the assumption we made for approach 1 for modeling the reaction. We can see that under the conditions we assumed, the reaction behaves like a Michaelis-Menten reversible reaction.

  3. Inhibition of directions:

    Examination of the in vitro properties of the recombinant AKR1C2 showed that it was potently inhibited in the oxidation direction by NADPH[6]. Since all AKRs catalyze an ordered bi-bi reaction, we can assume the enzyme we chose to use - AKR1C9 - has the same mechanism. In that case, if our model describes this mechanism accurately, this attribute should be reflected in it.

    Our model describes a general AKR enzyme. The difference between AKR1C2 and AKR1C9 is probably not in the expression of the rate reaction, but at the rate reaction's constants. Because these constants can vary in different environments even for the same strain of enzyme, it means that each direction of the reaction can be inhibited on different conditions (as in [5] for rat liver and rat skin). For that reason, we will need to make an assumption about the constants in order to show that the model describes it.

    We'll assume that \({K_1} < < {K_2}\) so that for \(\left[ {NADPH} \right] \approx \left[ {NADP} \right] > > {K_1}\),

    \(\frac{{\left[ {NADPH} \right]}}{{{K_1}}} > > \frac{{\left[ {NADP} \right]}}{{{K_2}}}\)

    From this assumption and equation (II) we will get:

    \[\begin{array}{l}\frac{{{V_{{m_f}}} \cdot \frac{{\left[ {NADPH} \right]}}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} - {V_{{m_r}}} \cdot \frac{{\left[ {NADP} \right]}}{{{K_2}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}}}}{{1 + \frac{{\left[ {NADPH} \right]}}{{{K_1}}} + \frac{{\left[ {NADPH} \right]}}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} + \frac{{\left[ {NADP} \right]}}{{{K_2}}} + \frac{{\left[ {NADP} \right]}}{{{K_2}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}}}} = \\ = \frac{{{V_{{m_f}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} - {V_{{m_r}}} \cdot \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}} \cdot \frac{{\left[ {NADP} \right]}}{{{K_2}}} \cdot \frac{{{K_1}}}{{\left[ {NADPH} \right]}}}}{{\frac{{{K_1}}}{{\left[ {NADPH} \right]}} + 1 + \frac{{\left[ {DHT} \right]}}{{{K_s}}} + \frac{{\left[ {NADP} \right]}}{{{K_2}}} \cdot \frac{{{K_1}}}{{\left[ {NADPH} \right]}} + \frac{{\left[ {3\alpha diol} \right]}}{{{K_P}}} \cdot \frac{{\left[ {NADP} \right]}}{{{K_2}}} \cdot \frac{{{K_1}}}{{\left[ {NADPH} \right]}}}}\mathop \approx \limits^{\begin{array}{*{20}{c}}{\left[ {NADPH} \right] > > {K_1}}\\{\frac{{\left[ {NADPH} \right]}}{{{K_1}}} > > \frac{{\left[ {NADP} \right]}}{{{K_2}}}}\end{array}} \\ \approx {V_{{m_f}}} \cdot \frac{{\frac{{\left[ {DHT} \right]}}{{{K_s}}}}}{{1 + \frac{{\left[ {DHT} \right]}}{{{K_s}}}}} = {V_{{m_f}}} \cdot \frac{{\left[ {DHT} \right]}}{{1 + \left[ {DHT} \right]}}\end{array}\]

    As we can see, our model can describe the inhibition of oxidation. With a simple assumption in regards to the relationship between two constants, we demonstrated that this reversible equation can behave like a single direction reaction.

Correlation with wet-lab results

Our team succeeded in over-expressing the 3\(\alpha \)-HSD enzyme in E. coli, and validated it's activity by measuring NADPH fluorescence. The results of the various measurements demonstrated certain behaviors which can be explained by our model. In this section we will simulate the model in order to see the correlation with our wet-lab results.

The rate equations for our system:

\(\begin{array}{l}\frac{{d\left[ {3\alpha - diol} \right]}}{{dt}} = \left[ {3\alpha - HSD} \right] \cdot \overbrace {\frac{{{K_{ca{t_f}}} \cdot \frac{{\left[ {NADPH} \right]}}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} - {K_{ca{t_r}}} \cdot \frac{{\left[ {NADP} \right]}}{{{K_2}}} \cdot \frac{{\left[ {3\alpha - diol} \right]}}{{{K_P}}}}}{{1 + \frac{{\left[ {NADPH} \right]}}{{{K_1}}} + \frac{{\left[ {NADPH} \right]}}{{{K_1}}} \cdot \frac{{\left[ {DHT} \right]}}{{{K_s}}} + \frac{{\left[ {NADP} \right]}}{{{K_2}}} + \frac{{\left[ {NADP} \right]}}{{{K_2}}} \cdot \frac{{\left[ {3\alpha - diol} \right]}}{{{K_P}}}}}}^{f\left( {\left[ {DHT} \right]\,\,,\,\left[ {NADPH} \right]\,,\,\left[ {3\alpha - diol} \right]\,,\,\left[ {NADP} \right]} \right)} - {\gamma _{{{\deg }_1}}} \cdot \left[ {3\alpha - diol} \right]\\\frac{{d\left[ {NADP} \right]}}{{dt}} = f - {\gamma _{{{\deg }_2}}} \cdot \left[ {NADP} \right]\\\frac{{d\left[ {DHT} \right]}}{{dt}} = - f - {\gamma _{{{\deg }_3}}} \cdot \left[ {DHT} \right]\\\frac{{d\left[ {NADPH} \right]}}{{dt}} = - f - {\gamma _{{{\deg }_4}}} \cdot \left[ {NADPH} \right]\end{array}\)

We will solve the equation system numerically using Matlab with the constants from the table below:

Parameter Value Units source comment
\({{K_{ca{t_f}}}}\) \({\rm{0}}{\rm{.563}} \cdot {\rm{1}}{{\rm{0}}^{ - 3}}\) \(\frac{1}{{\sec }}\) Calculated from [5] using the molar mass of AKR1C9 from [3] (***)
\({{K_{ca{t_r}}}}\) \({\rm{1}}{\rm{.628}} \cdot {10^{ - 3}}\) \(\frac{1}{{\sec }}\) Calculated from [5] using the molar mass of AKR1C9 from [3] (***)
\(\left[ {3\alpha - HSD} \right]\) 0.1 \(\mu M\) Chosen arbitrarily. The concentration of the enzyme is connected linearly to the maximum reaction rates and do not effect the dynamics.
\({K_s}\) 0.38 \(\mu M\) From [5] (***)
\({K_p}\) 2.79 \(\mu M\) From [5] (***)
\({K_1}\) 0.38 \(\mu M\) None (****)
\({K_2}\) 2.79 \(\mu M\) None (****)
\({\gamma _{{{\deg }_1}}}\) 0 \(\frac{1}{{\sec }}\) None (*)
\({\gamma _{{{\deg }_2}}}\) 0.0001 \(\frac{1}{{\sec }}\) None (**)
\({\gamma _{{{\deg }_3}}}\) 0 \(\frac{1}{{\sec }}\) None (*)
\({\gamma _{{{\deg }_4}}}\) 0.0001 \(\frac{1}{{\sec }}\) None (**)

(*) Note that since we could not measure DHT and \(3\alpha-diol\), we neglected them in the simulation for the sake of simplicity.

(**) The degradation rates we measured for NADPH varied too much between measurements, probably due to the different concentration of lysates. Since the purpose of this model is to emulate the behavior of our system, small degradation constants were chosen arbitrarily for NADPH and NADP. By giving both cofactors small degradation constants, we are able to mimic the behavior of the reaction while looking at the dynamics of the enzyme's reaction.

(***) We chose to use the constants measured at [5], although they are probably not the correct constants for the environment of the experiment. For that reason, these constants sometime differ from simulation to simulation by orders of magnitude in order to get the simulation result that best fit our system's behavior.

(****) Since the values of \({{K_1}}\) and \({{K_2}}\) are unknown, and since we demonstrated how the relation between them can effect the reaction rate (inhibition of direction), we chose an arbitrary values for most simulations that are identical to \({{K_s}}\) and \({{K_P}}\). That being said, we experimented a lot with these values in order to test different behaviors of this system.


  1. Activity check simulation

    We simulated the system in time domain under the initial conditions:

    \[\left( {\left[ {DHT} \right],\,\left[ {NADPH} \right],\,\left[ {3\alpha - diol} \right],\,\left[ {NADP} \right]} \right) = \left( {{{40}_{\left[ {\mu M} \right]}},\,{{150}_{\left[ {\mu M} \right]}},\,0,\,0} \right)\]

    which are the same initial conditions for the avctivity check experiment

    Figure 5 – NADPH concentration vs. Time. We can see that the behavior of this simulation is close to the result of the activity check.

  2. Reaction rate vs. DHT concentrations

    We simulated the system for varying initial concentrations of DHT at \[\left[ {NADPH} \right] = {150_{\left[ {\mu M} \right]}}\]. For each initial concentration we took the derivative of NADPH at t=0 (in absolute value), which is the maximal reaction rate for the set of initial conditions.

    Figure 6 - Initial reaction rate vs. Initial DHT concentration. The behaviour is similar to the results.

    The result from the simulation is very similar to the wetlab result, with the reaction rate saturating as the concentration of DHT increases. This can be explained by our model since we saw that for low concentrations of \(3\alpha - diol\) and NADP the reaction is close to michaelis menten.

  3. Enzymatic activity as a function of NADPH concentration

    As noted in the results section, we observed a decrease in the reaction rate with an increase in NADPH concentration in presence of 3ɑ-HSD enzyme, contrary to substrate dependency. One possible explanation for this is cofactor inhibition - the addition of DHT to the lysate with the enzyme and NADPH was done 30 minutes after adding the NADPH to the lysate, in order for the reaction between NADPH and the lysate to reach steady state. During that time, there was a decrease in NADPH levels from their initial values. We hypothesized that some of the NADPH was converted to NADP, which in turn was bound to the enzyme. Since There is no \(3\alpha-diol\) in the system, the oxidation reaction cannot take place so the enzymes that are bound to NADP molecules cannot participate in another reaction. In this state, NADP inhibit the reduction of DHT.

    In order to verify that the model can actually describe this result, we took the initial concentrations of NADPH from the experiment and checked those levels under 30 minutes. We made an assumption that most of the molecules of NADPH that broke down before the addition of DHT were converted to NADP. Under this assumption, we ran a simulation with the initial conditions of the wetlab experiment after 30 minutes:

    Figure 7 - Initial reaction rate vs. Initial NADPH and NADP concentrations.

    As can be seen in the figure, the derivatives were the high for the reactions in which the NADPH/NADP ratio was high, and vice versa.

    In order to better understand the process and the results, a simulation of the reaction rate of NADPH as a function of both initial NADPH and initial NADP was built. In each run we changed the ratio between K1 and K2. The result is in the following video:

    from this video we can see that for K1>K2, the system is more sensitive to an increase in NADP (which inhibit the reaction) than to an increase in NADPH.

    This can explain the results we got in the lab because the more we increased the initial concentration of NADPH, the more NADPH was converted to NADP and so the inhibition effect was more significant.


1. Cooper, W. C., Heredia, V. V., Jin, Y., & Penning, T. M. (2006). Comparison of the Rate-Limiting Steps in 3alpha-Hydroxysteroid Dehydrogenase (AKR1C9) Catalyzed Reactions. In H. Weiner, B. Plapp, R. Lindahl, & E. Maser, Enzymology and Molecular Biology of Carbonyl Metabolism 12 (pp. 301-307). Purdue University Press. Retrieved from
2. Jin, Y., Heredia, V. V., & Penning, T. M. (2006). Enzymatic Mechanism of 5alpha-Dihydrotestosterone Reduction Catalyzed by Human Type 3 3alpha-Hydroxysteroid Dehydrogenase (AKR1C2): Molecular Docking and Kinetic Studies. In H. Weiner, B. Plapp, R. Lindahl, & E. Maser, Enzymology and Molecular Biology of Carbonyl Metabolism 12 (pp. 294-300). Purdue University Press. Retrieved from
3. Penning TM, M. I. (1984, Sep 15). Purification and properties of a 3 alpha-hydroxysteroid dehydrogenase of rat liver cytosol and its inhibition by anti-inflammatory drugs. The Biochemical Journal. Retrieved from
4. Penning, T. M., Jin, Y., Heredia, V. V., & Lewis, M. (2003). Structure–function relationships in 3alpha-hydroxysteroid dehydrogenases: a comparison of the rat and human isoforms. The Journal of Steroid Biochemistry & Molecular Biology, 85, 247-255. Retrieved from
5. Pirog, E. C., & Collins, D. C. (1994, April). 3alpha-Hydroxysteroid dehydrogenase activity in rat liver and skin. Steroids, 59, 259-264. Retrieved from
6. TEA LANISˇ NIK RIZˇ NER*, HSUEH K. LIN, DONNA M. PEEHL, STEPHAN STECKELBROECK,DAVID R. BAUMAN, AND TREVOR M. PENNING. Human Type 3 3alpha-Hydroxysteroid Dehydrogenase Aldo-Keto Reductase 1C2) and Androgen Metabolism in Prostate Cells. Retrieved from
7. Biswas, M. G., & Russell, D. W. (1997). Expression cloning and characterization of oxidative 17β-and 3α-hydroxysteroid dehydrogenases from rat and human prostate. Journal of Biological Chemistry, 272(25), 15959-15966. Chicago
8. Rob Phillips Jane Kondev (Author, Julie Theriot, (2008) "Physical Biology of the Cell" p.219-237, 978-0815341635

Contact Us