Modelling our system:

For a better development of our system we thought of modelling. The possibility of concentration measurement with our system needs a way of calibration. Therefore the antigen concentration in the sample needs to correlate with the color shift ratio in such a way, that a calibration line in the desired concentration range can be determined. This could be predefined by the ratio of the three compounds. E.g. having more detection component could lead to a higher sensitivity for higher antigen concentrations but could also be less sensitive to very low antigen amounts, since the background noise could be higher and the antigen concentration is being overlooked by this.

The possible quantification should be based on a calibration based on the correlation between the time and the antigen amount. Therefore we define a threshold value of the colour-shift ratio and measure/simulate the time that is needed for a sample (with a specific antigen amount) to reach this threshold. The more this times can be discriminated, the better.

For the modelling of the simulation we need to obtain time dependent values of the color ratio based on some constants (to be determined) and the initial amount of the components (can be adjusted for every simulation run). Since the protease is an enzyme, we use the Michaelis-Menten-kinetic to modell our system. For our first model engineering cycle we simply had the formula:

d[P]/dt = v = [S]/([S]+KM) * vmax.

explanation: [P] = product concentration, d[P]/dt = change of Product concentration over the change of the dt amount of time, v = enzymatic rate, [S] = substrate concentration, KM = Michaelis-Menten-constant, vmax = maximal enzymatic rate

We thought thereby of [P] as the concentration of our protease-cutted reporter component and [S] as the concentration of our protease-uncutted Reporter- and Amplification component. But this is an too sketchy model since the Reporter component is not only contributing to [S] (since it has an protease recognition site) but also is by itself having an protease activity once activated and thereby contributing to a raising vmax over the time. To be able to model this we can have a further look at vmax. It can be written as: vmax = kcat * [E0]. Thereby [E0] is the totally available enzyme concentration. For our second model engineering cycle we therefore assume, that our amplification component is contributing to the [E0] and the [S]. Thats why we came up to split the formula into the following two formulas:

I: d[P]/dt = v_rep = [S_rep]/([S_rep]+KM) * kcat * [E0]

explanation: [P] = concentration of cutted reporter component, v_rep = enzymatic rate of cutting the reporter component, [S_rep] = uncutted reporter component concentration

II: d[E0_amp]/dt = v_amp = [S_amp]/([S_amp]+KM) * kcat * [E0]

explanation: [E0_amp] = concentration of the cutted amplification component, v_amp = enzymatic rate of cutting the amplification component, [S_amp] = uncutted amplification component concentration

We now have the amplification and reporter component in our model. But what is with the detection component? We made the assumption, that together with the background noise we can set it as the initial [E0]. Usually we believe the behaviour of the detection component to be in a limited growth fashion dependent on the actual concentration of the detection component and the amount of antigen present. But to make the model not to complex we simply set this value at a constant level. Additionally we wanted to make a proof of concept of our system using only the TEV protease replacing the detection component. Therefore we thought an initial [E0] is a good option. Following assumption is bringing in the Detection component into the context of the other formulas:

[E0] = [E0_det] + [E0_amp]

explanation: [E0_det] = concentration of the active detection component

For our third model engineering cycle we brought in the time since we want to have a look of the colour shift over time. Therefore we thought of a numeric approach. By using small time difference values for dt we can approximate the behaviour of our numbers over time. This is also known as the Gaussian method named after the famous mathematician who originated in Braunschweig, as we do too :D

Therefore we needed a third formula explaining the change of Substrate over time. This third formula is based on the assumption, that the product formation is stoichiometric equivalent to the substrate depletion:

III: d[S]/dt = [S_rep]/dt + [S_amp]/dt = - d[P]/dt - d[E0_amp]/dt

Know we can define the start values for t0 and for t+dt, the dt dependent change of the values:

t0: (explanation: the initial defined values)

[P](t=0) = 0 (explanation: at t=0 no product was fromed so far)

[E0](t=0) = [E0_det]

[S](t=0) = [S0] = [S0_rep] + [S0_amp] (explanation: the initial substrate concentration for the reporter and the amplification component)

t+dt:

[P](t+dt) = [P](t) + v_rep(t) * dt (solved with the help of formula I)

[E0](t+dt) = [E0](t) + v_amp(t) * dt (solved with the help of formula II)

[S](t+dt) = [S](t) - d[P](t) - d[E0](t) (solved with the help of formula III and the simplifications if I and II: d[P](t) = v_rep * dt as well as: d[E0](t) = v_amp * dt)

This results to following graph:

But this looked not very realistic. Especially when the maximal product concentration is reached whe would find it more realistic if the buckling would be more curved. So we made a new model based on the Michaelis-Menten kinetic. this time with the enzyme substrate complex. We also put a new kind of Substrate within our formulas: A. The product of A will be enzyme E. This is representing our amplification component. We made the following differential equations the basis for our new model:

d[E]/dt = - kf [E] [S] + kr [ES] + kcat [ES] - kfa [E] [A] + kra [EA] + 2 kcata [EA]

d[S]/dt = - kf [E] [S] + kr [ES]

d[ES]/dt = kf [E] [S] - kr [ES] - kcat [ES]

d[P]/dt = kcat [ES]

d[A]/dt = -kfa [E] [A] + kra [EA]

d[EA]/dt = kfa [E] [A] - kra [EA] - kcat [EA]

Following curves were simulated:

This sigmoidal curve looks more realistic to us. unfortunately now we have many additional constants that are difficault to obtain from characterising a part. Nevertheless we made another engineering cycle. We wanted to bring in the time that is needed for the cutted amplification component to get active (i.e. the time needed for the cage to dissociate, the split intein halves to reconstitute, the reconstituted intein to splice the TEV protease and the reconstituted TEV protease to became active). Therefore we build in another compound, the first product of the enzyme-amplification-component complex. This is in a single way reaction further processed to the TEV protease. The turn over point of this reaction is thought to be slow since it will represent the time that is needed for the cutted amplification component to have an active TEV activity. Following equation was changed compared to the beforehand formulas:

d[E]/dt = - kf [E] [S] + kr [ES] + kcat [ES] - kfa [E] [A] + kra [EA] + kcata [EA] + kcatpa [PA]

An Additional formula is needed, describing the concentration of the intermediate:

d[PA]/dt = kcata [EA] - kcatpa [E]

With this new model we gained even more realistic curves:

We can see, that the first exponential phase is slower then in the previous model since some TEV protease activity needs first to be gained by activation of the amplification component. But the slowered product formation near the saturation changes not as slow as the rising product formation rate at the beginning. This is an awaited behaviour based on the mode of action of the amplification component within our system.

With this optimised model we can now compare the influence of altered compound compositions on the effect of our simulated results. Therefore we compare two simulations:

The third simulation was conducted with a concentration of the amplification component of [A0] = 2 mol/l while in the fourth simulation a concentration of [A0] = 5 mol/l was used. For both simulations the antigen concentrated was simulated with [E0] = 0.01 and [E0] = 0.001 for the high and low antigen concentration respectively.

The threshold value was set to 4. So the time after [P] >= 4 was reached was extracted and the time difference between the high and the low antigen amount was calculated. A higher difference can lead to a better calibration line. We got a difference of 25.6 seconds for the third and 26.3 seconds for the fourth simulation. Therefore we can suggest to use more ampification compound to gain a better calibration line! To see how we exactly made a simulation and how we generated the graphs. Please have a look further down on this page, at the appendix section. There you can find the python code.

Outlook:

In this model we definitly need more realistic values. Therefore the parts should be characterised further. But with the establishement of this model we could possible find out the best compositions of our compounds for our system to work at its best. This can save much time consuming and blunt lab work.

Appendix:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created 2022

@author: Samuel Nestor Meckoni
"""

import numpy as np

import matplotlib.pyplot as plt


# concentration values for t = 0
E_0 = 0.01
P_0 = 0.0
S_0 = 10
A_0 = 2

ES_0 = 0
EA_0 = 0
PA_0 = 0

# time scaling
t = np.linspace(0,100,1000)
dt = t[1] - t[0]

# constants
kcat = 10
kf = 0.1
kr = 0.1
kcata = 9
kfa = 0.09
kra = 0.09
kcatpa = 0.1


def iteration( E_0, P_0, S_0, A_0, ES_0, EA_0, PA_0, t, dt, kcat, kf, kr, kcata, kfa, kra, kcatpa ):
    # calculates the concentration changes over time and results tables
    
    for i in range(len(t)):
    #for loop for each iteration
        if t[i] == 0:
        # initialising values are written in tables/arrays for t=0
            E_t = [E_0]
            P_t = [P_0]
            S_t = [S_0]
            A_t = [A_0]
            ES_t = [ES_0]
            EA_t = [EA_0]
            PA_t = [PA_0]
        else:
        # Here are the actual values calculated, corresponding to the formulars.
            E_t.append( E_t[i-1] + ( dt * ( -kf * E_t[i-1] * S_t[i-1] + kr * ES_t[i-1] + kcat * ES_t[i-1] -kfa * E_t[i-1] * A_t[i-1] + kra * EA_t[i-1] + kcata * EA_t[i-1] + kcatpa * E_t[i-1] ) ) )
            P_t.append( P_t[i-1] + ( dt * ( kcat * ES_t[i-1] ) ) )
            S_t.append( S_t[i-1] + ( dt * ( -kf * E_t[i-1] * S_t[i-1] + kr * ES_t[i-1] ) ) )
            A_t.append( A_t[i-1] + ( dt * ( -kfa * E_t[i-1] * A_t[i-1] + kra * EA_t[i-1] ) ) )
            ES_t.append( ES_t[i-1] + ( dt * ( kf * E_t[i-1] * S_t[i-1] - kr * ES_t[i-1] - kcat * ES_t[i-1] ) ) )
            EA_t.append( EA_t[i-1] + ( dt * ( kfa * E_t[i-1] * A_t[i-1] - kra * EA_t[i-1] - kcata * EA_t[i-1] ) ) )
            PA_t.append( PA_t[i-1] + ( dt * ( kcata * EA_t[i-1] - kcatpa * E_t[i-1] ) ) )
    return E_t, P_t, S_t, A_t, ES_t, EA_t, PA_t
          

# call of the function  
E_t, P_t, S_t, A_t, ES_t, EA_t, PA_t = iteration( E_0, P_0, S_0, A_0, ES_0, EA_0, PA_0, t, dt, kcat, kf, kr, kcata, kfa, kra, kcatpa )
r = P_t
E_0 = 0.001
E_t, P_t, S_t, A_t, ES_t, EA_t, PA_t = iteration( E_0, P_0, S_0, A_0, ES_0, EA_0, PA_0, t, dt, kcat, kf, kr, kcata, kfa, kra, kcatpa )
s = P_t

# setting new values for the next two simulations
A_0 = 5
E_0 = 0.01
E_t, P_t, S_t, A_t, ES_t, EA_t, PA_t = iteration( E_0, P_0, S_0, A_0, ES_0, EA_0, PA_0, t, dt, kcat, kf, kr, kcata, kfa, kra, kcatpa )
r2 = P_t
E_0 = 0.001
E_t, P_t, S_t, A_t, ES_t, EA_t, PA_t = iteration( E_0, P_0, S_0, A_0, ES_0, EA_0, PA_0, t, dt, kcat, kf, kr, kcata, kfa, kra, kcatpa )
s2 = P_t


# creation of the figures
fig, ax = plt.subplots()
ax.plot(t, s, color='b', label='low antigen amount')
ax.plot(t, r, color='r', label='high antigen amount')
legend = ax.legend(loc='upper left', shadow=False, fontsize='x-large')
legend.get_frame().set_facecolor('C6')
ax.set(xlabel='time in s', ylabel='[P] in mol/l',
       title='Third Simulation run.')
ax.grid()
fig.savefig("graph4.png")
plt.show()

fig2, ax2 = plt.subplots()
ax2.plot(t, s2, color='b', label='low antigen amount')
ax2.plot(t, r2, color='r', label='high antigen amount')
legend = ax2.legend(loc='best', shadow=False, fontsize='x-large')
legend.get_frame().set_facecolor('C6')
ax2.set(xlabel='time in s', ylabel='[P] in mol/l',
       title='Fourth Simulation run.')
ax2.grid()
fig2.savefig("graph5.png")
plt.show()


# calculation of the threshold value (i.e. the time, at when a [P] of 4 was first reached)
rsr2s2=[]
for i in range(len(t)):
    if r[i] >= 4:
        rsr2s2.append(t[i])
rsr2s2 = [rsr2s2[0]]
for i in range(len(t)):
    if s[i] >= 4:
        rsr2s2.append(t[i])
rsr2s2 = rsr2s2[0:2]
for i in range(len(t)):
    if r2[i] >= 4:
        rsr2s2.append(t[i]) 
rsr2s2 = rsr2s2[0:3]
for i in range(len(t)):
    if s2[i] >= 4:
        rsr2s2.append(t[i])
rsr2s2 = rsr2s2[0:4]
print(rsr2s2)
print("Time difference for the third simulation:",rsr2s2[1]-rsr2s2[0])
print("Time difference for the fourth simulation:",rsr2s2[3]-rsr2s2[2])