Team:SYSU-CHINA/Try







Modeling


How many Dox, How cell response and others problems
How many Dox, How cell response and others problems

Our project aims to regulate the activity of CAR-T cells in order to reduce the risk of inflammatory factor storm which do damage to patients. In order to further reveal how many drugs and how often we should give patients, we have modeled to solve these problems. In order to guide the amount of administration, we expect to be able to know how the amount of administration or the time after administration affect the effect of drug, in order to avoid the calculation trouble caused by the excessive definition of the variable, we divided the whole process into ①dox pharmacokinetics of Dox②Tet-On system‘s response and U24 protein expression ③how U24 protein affect concentration of CAR and ④the impact of cell surface concentration of CAR and ligand on the immune intensity.

1、Why we use dox

Dox is a tetracycline analogue that induces Tet-On system activation. Compared with tetracycline, it has better stability and lower toxicity. Therefore, we choose Dox as our inducer.

2、Why we choose CAR as the points of regulation?

Observing the structure of the entire T cell activation signaling pathway and selecting the DNA-RNAP complex for sensitivity analysis, it can be easily seen that only the sensitivity to CAR is composed of positive terms only, so the sensitivity to CAR should be the highest.


Figure1:pathway of T-Cell activation[1]
3、 Concentration in cells

As we all known，drugs are not all distributed in cells but are distributed between plasma and cells in human‘s body’. We consulted the literature and obtained the transfer rate constant of drugs form the digestive system to plasma or form cytosol to the plasma as well as redistribution rate constant form cytosol to the plasma [2]，[3]. After consider the drug’s decomposition rate is consistent in the plasma and cell, we list reaction equation as follows


By solving the homogeneous ODEs（ordinary differential equations ）obtained by the law of mass action, we finally get the function of the concentration of the drug in the cytosol about time (see in: supplementary information) and plot it as follows(assume dosage is 2mmol/kg). As can be seen from the analysis of the Dox -time relationship, Dox concentration-Time relationship follows a similar manner to a three-exponential function will reach a peak after about 8 hours, and then decrease quickly, the concentration of dox has dropped very low at about 4000 minutes. We also plot the situation of successive administration and find the drug concentration in the cytosol will stabilize when administered 2mmol per kg per 1000 minutes with certain parameters we use.


Figure2: time dependent cytosol Dox concentration



Figure3: situation of successive administration 2mmol Dox per kg per 1000 minutes
4、Response of Tet on system

The change of dox over time is obtained. Next step is to consider the Response of the Tet-on system. Due to the complexity of the above functions, we can not find the analytical solution so we assume that under normal conditions, the content of rtTA in the cells is in steady state and the transcription rate is determined by the following formula:.
v[t]=B⁄((1+dox[t]/kd))+A (((dox[t])⁄kd))⁄((1+dox[t]/kd))
Among them, (1+dox[t]/kd) and (((dox[t])⁄kd))⁄((1+dox[t]/kd)) represent the ratio of time-dependent Dox-unbound /binding DNA, B and A are the rate of leaked transcription and the rate of induced transcription, respectively, and kd is the apparent dissociation constant. We implicitly include rtTA and some other constant in kinetic constants above according to literature and experiments. Having obtained certain data [4]、[5]、[6]、[7]., we can get certain transcription rate and solve the differential equation with considering the degradation of RNA which is listed below.
rna'[t]=v[t]-kdeRNArna[t],rna[0]=0


Figure4 : RNA concentration-Time relationship

After considering the translation rate constant by estimated from the kinetic constant and the ribosome concentration[7],[8],[9], we try to solve the ODE about translation. However, due to the use of interpolation function, we can't solve the equation. At first, we want to solve it by the perturbation method which mean determine the protein yield without degradation firstly and use the result to estimate the degradation rate at each moment, but we found that it is difficult to construct a suitable perturbation function over a long period of time because the value of the perturbation always become too colossal, thereby we having to give up this method.
Next, we remember that we can use the exponential decay law to estimate this result, so we estimate the yield-time relationship according to the expression, and find the result qualitatively.
u24[t]=E^(-kdeu24 t )Integrate[ktrans* rna[t],t]
Qualitatively, the model is correct, but the drawing of the rate of production shows that such treatment overestimates the degradation. We use the average time, which is obtained by integrate the amount of protein without considering degradation with time and divide the result by the amount of protein, instead of the time to improve calculation result but it has not improved much. We give the relationship between degradation rate constant and u24 concentration-Time relationship, too. We find that the degradation rate constant affect maximum U24 concentration and retention time in great extent. It may possible to use variant degradation tag to control concentration of U24 and density of CAR in wide range.


Figure5: u24 concentration-Time relationship



Figure6: u24 concentration-Time relationship with consider degradation dependent on average time



(a)
(b)

(c)
(d)
Figure7: Effect of degradation rate constant on and U24 Concentration-Time relationship
(a)kde=0 (a)kde=0.003(a)kde=0.012(a)kde=0.0332

After surveying the equation and using the new micro-element ktrans rna[t] E^(-kdeu24(t0-t)dt instead of dy and integrate, we get an integral upper bound function. After finding the numerical solution of it and plot it as follows, we find those attenuation result is quite accuracy, RNA degrade slowly and lasts longer. Further calculation will be obtained after we get higher computing power. In general, the subsequent refinement changes quantitatively, but has little effect on the peak shape, and is in good agreement with the tendency of amount of dox.


F igure8: u24 concentration-Time with change integrate micro-element
5、Hill equation

Because we have some doubts about the accuracy of some constants use above, we hope to use Hill equation to get another new model. After obtain the protein yield-time relationship by experiment, we find that it roughly matches the Hill equation. In the form, we use the Hill equation too fit it which is often used. Since cytoplasm Dox concentration can be regarded as a constant during cell culture, according to the experimental data, expression
exp[t_]:=A1 (dox[t]^ h)⁄((dox[t]^h+ki^h))
is fitted to obtain a function of u24 over time, with fewer guess components than we obtain in above chapter. In this expression, A1 stand for the production of induced express, ki stand for affinity constant between dox and DNA, h is Hill constant. Data is obtained by our experiment about cell-culture with concentration gradient Dox.


Figure9: u24 concentration-Time determined by Hill equation and experiment

The expression is
1.13701*(Dox^0.924472)⁄((Dox^0.924472+0.0731035921466379^0.924472))
Obtained the Hill equation about Dox and u24, we use Dox[t] which is obtained above instead constant Dox concentration and get a new graph about U24 and time.

6、the relationship between car and immune signal intensity
Why choose A3 as an indicator?

We know that in immunology, surface receptors are generally transmit information to the downstream pathway through phosphorylation. Due to the phosphorylated region on the TCR is mainly the ITAM domain and there are three ITAM domains on the CAR, it is reasonable to use the concentration of CAR which three ITAMs are phosphorylated as the indicator of immune signal intensity.

Dynamic mode diagram

We chose a model called “separation model” to describe the behavior of CAR during the immunization process. The main content of this model is that large phosphatases will be "squeezed" out of the immune synapse to form a kinase-rich domain(KRD) during the process of immunization process, T cell receptors that bind to pMHCs will undergo a series of phosphorylation and initiate downstream pathways by recruiting molecules such as ZAP70. After reference to the dynamic proofreading model, we made a dynamic mode diagram of CAR phosphorylation.[10],[11],[12],[13],[14]
According to the law of mass action, we list a series of ordinary differential equations(ODE) as follows to describe phosphorylation and use the density of the three-phosphorylated receptor(A3) as an indicator of immune intensity. Since the ordinary differential equations have no analytical solution, we adopt two approaches to solve this problem.


ODEs:
T^' [t]=-konT[t]*M[t]+koff(a1[t]+a2[t]+a3[t]+a0[t])
M^' [t]=-konT[t]*M[t]+koff(a1[t]+a2[t]+a3[t]+a0[t])
a0^' [t]=konT[t]*M[t]-koffa0[t]-pa0[t]+qa1[t]
a1^' [t]=pa0[t]-qa1[t]+qa2[t]-pa1[t]-koffa1[t]
a2^' [t]=pa1[t]-qa2[t]+negaa3[t]-posia2[t]-koffa2[t]
a3^' [t]=posia2[t]-negaa3[t]-koffa3[t]
T[0]=100,M[0]=300,a0[0]=0,a1[0]=0,a2[0]=0,a3[0]=0
Numerical Solution

Using matlab, we get the numerical solution of the system of ODEs. It can be seen that the tertiary phosphorylated receptors gradually increased and eventually almost all of the receptors were phosphorylated. We also made a diagram about how the cell surface CAR density affects the density of the phosphorylated receptors.


Figure10：A3-t



Figure11：T0-M0-A3
Steady state equation:

We assume that all phosphorylated receptors are in steady state during the immune response which changes ODEs into a homogeneous equation, making it easier to solve. Using Mathematica, we express A3 as a function of CAR receptor density T0 and its ligand density M0 (see in supporting information) and draw a three-dimensional graph of it. It can be seen that A3 is almost linearly positively correlated for T0 and M0, respectively. Taking the ratio of the density of A3 produced as a specificity index, we plot how T0 density affects CAR-T’s specificity of cells have different ligand numbers and see when ratio of density of ligands between two cells is low, specificity depends on density of T0 linearly unless T0 lower than a certain number, which will make specificity vanish. When the ratio approach or exceed 0.3, specificity vanish totally. Using mathmatica, we also plot how kinetic parameter affect our result.


Figure12：Effects of CAR receptor density T0 and its ligand density M0 on A3



(a)
(b)

(c)
(d)
Figure13：effects of cell surface T0 density on specificity of ligand density（r stand for ratio of density of ligands between two cells(a)ratio between decoy cell and target cell equal to 0.004(b)ratio between decoy cell and target cell equal to 0.094(c)ratio between decoy cell and target cell equal to 0.232(d)ratio between decoy cell and target cell equal to 0.306


Figure14: effects of kinetic parameters on A3(T0，M0)，kinetic parameters are show above
Monte Carlo and diffusion

The diffusion of receptor between the kinase-rich domain and the phosphatase-rich environment is important. We use the Monte Carlo algorithm and the random walking model to simulate the diffusion of the CAR receptor on the cell surface and calculate the average retention time and equilibrium density ratio of those two regions. After add a diffusion-related term to the ODEs and chemical equation, we get a new Numerical Solution and find an interesting result——there is no steady state in this model, instead, production of A3m show a peak shape after the formation of immune synapses, and then returns to 0. We also used this model to do some simulations of KRD formation. New chemical equations we use are show below.






Figure15: Result of simulation diffusion by Monte Carlo algorithm and the random walking mode
Figure16：A3-t relationship with consider diffusion
Figure17：a short movie simulates CAR on cell surface during T-cell activation，in which blue point stand for A0 and red point stand for A3

According to the literature, phosphorylated CAR will form cluster on the cell surface, and then zap70 will be recruited to initiate the downstream reaction. Here, we use ZAP70 concentration to more accurately reflect the immune intensity. To make calculate easier, we define the biggest cluster as a tetramer and assuming that each CAR recruits only one ZAP70, according to the aggregation- recruitment pattern and the low of mass action, and using the steady-state approximation, we solve the expressions of the different state complexes and determine the relationship of the amount of activated ZAP70. T0, M0 with certain various kinetic constants, the final result is too complicated to listed here[10]


Figure18: schematic diagram of more accurately model(A) Schematic of the model(B) State space of an individual cluster. Each state is identified by the doublet of integers ðn;mÞ: the number of ligated receptors (n) and the number of those receptors that are ZAP70 positive (m). Note that m<=n.(C) Possible stochastic transitions, and associated rates, from a representative state (n;m)[15]
7、Conclusion

Based on the relationship between dox and A3, we can get a new function for the dose and actual effect and understand the relationship between dose, administration time and immune intensity better. Combined with more clinical trials, this new function will be more instructive. In combination with relationship between Dox and immune intensity, we can get guidance on the dose of dox to the patient, too.
Using the modeling we modeled, we can learn that Dox should be 8 hours administered before the symptoms worsen. If long-acting inhibition is required, it should be administered at a frequency of about 1 time/day (based on the dose of 5 mmol/kg). The combination of models and more experimental studies in the clinic, we can better guide the application.

Reference

[1] KEGG, PATHWAY: T cell receptor signaling pathway map04660[2018.10.03] https://www.genome.jp/dbget-bin/www_bget?map04660
[2] Pharmacokinetics of doxycycline hydrochloride in cherry valley duckling during aflatoxicosis , Xiaoli He, Gang Ye ,et.al., Int. J. Clin. Exp. Med 2016;9(7):13035-13040
[3] Pharmacokinetics and bioavailability of doxycycline hyclate after oral administration in calves, L A. Meijer , K. G F. Ceyssens, et.al., Veterinary Quarterly, 1993; 15: 1-5
[4] Construction and Modelling of an Inducible Positive Feedback Loop Stably Integrated in a Mammalian Cell-Line, Velia Siciliano, Filippo Menolascina, et.al., PLoS Computational Biology, June 2011 | Volume 7 | Issue 6 | e1002074
[5] Using the Tet-On system to develop a procedure for extracting transcription factor activation dynamics. Z. Huang, C. Moya, A. Jayaraman, and J. Hahn. Molecular bioSystems, 6(10):1883{1889, October 2010.
[6] Stochastic simulations of the tetracycline operon, Konstantinos Biliouris, Prodromos Daoutidis, Yiannis N Kaznessis, BMC Systems Biology 2011, 5:9
[7]: Stochastic Kinetic Analysis of Developmental Pathway Bifurcation in Phage l-Infected Escherichia coli Cells, Adam Arkin, John Ross and Harley H. McAdams, Genetics 149: 1633–1648 (August 1998)
[8]: Kinetic dissection of fundamental processes of eukaryotic translation initiation in vitro, Jon R.Lorsch1,2 and Daniel Herschlag2, The EMBO Journal Vol.18 No.23 pp.6705–6717, 1999
[9]: Kinetic and thermodynamic analysis of the role of start codon/anticodon base pairing during eukaryotic translation initiation, SARAH E. KOLITZ, JULIE E,et.al. RNA (2009), 15:138–152
[10] Ligand Detection and Discrimination by Spatial Relocalization: A Kinase-Phosphatase Segregation Model of TCR Activation , Nigel J. Burroughs, Zorana Lazic, and P. Anton vander Merwey, Biophysical Journal Volume 91 September 2006 1619–1629
[11] SEGREGATION MODELS, Elaine P. Dopfer, Mahima Swamy,et.al. Multi-chain Immune Recognition Receptor Signaling: From Spatiotemporal Organization to Human Disease, Springer Science & Business Media, Dec 16, 2008
[12] KINETIC PROOFREADING MODEL Byron Goldstein, Daniel Coombs.et.al, Multi-chain Immune Recognition Receptor Signaling: From Spatiotemporal Organization to Human Disease, Springer Science & Business Media, Dec 16, 2008
[13] SERIAL TRIGGERING MODEL Jacob Rachmilewitz,,Multi-chain Immune Recognition Receptor Signaling: From Spatiotemporal Organization to Human Disease, Springer Science & Business Media, Dec 16, 2008
[14] Modelling T Cell Activation, Cliburn Chan Chi Wei, Centre for Nonlinear Dynamics and its Applications University College London, January 14, 2002 [P]
[15] A DNA-Based T Cell Receptor Reveals a Role for Receptor Clustering in Ligand Discrimination, Marcus J. Taylor, Kabir Husain,et.al. Cell 169, 108–119,March 23, 2017