Team:UCAS-China/Model
MODEL
Overview
Our model consists of three parts.
In the first part, we establish OED-based models to verify the superiority of our design of cold-induced gene circuit, and to describe and predict the function of our gene circuit.
In the second part, we employ an OED-based model to compare the distribution of drug by two different ways of drug delivery--pills and bacteria, and illustrates that delivering drug by our bacteria is better than traditional method.
The third part is about our hardware. We simulate the temperature distribution in intestinal tract when the capsule heats. We also implement the location algorithm of magnetic capsule based on Levenberg-Marquardt Method.
Part 1 Cold-induced Thermosensitive System
In this part we establish ODE-based models for the cold-induced thermosensitive system. Firstly, we establish model for TEV-TF subsystem, and assist in designing of gene circuit by valuing two different designs. We then establish a model for mflon-GOI subsystem and assess the degradation rate. We finally determine some parameters to describe the system by fitting the data from the experiment more precisely
1.TEV-TF subsystem
Parameter used in this model：
| $$g_{x}$$ | Generation rate of x |
| $$d_{x}$$ | degradation rate of x |
| $$C_{x}$$ | Concentration rate of x |
| $$\alpha_{x}$$ | Self-degradation rate of x |
| $$\beta_{x}$$ | Thermosensitive coefficient of x |
| $$n_i$$ | Hill coefficient |
| $$K_h$$ | K for hill equation |
| $$K_m$$ | K for Michaelis-Menten equation |
| $$K_h$$ | K for hill equation |
| $$V_{TEV,max}$$ | Max generation rate of TEV |
| $$V_{TF,max}$$ | Max degradation rate of TF |
Figure 1. Genetic circuit of TEV-TF subsystem
For TEV:
$$\frac{d C_{T E V}}{d t}=g_{T E V}-d_{T E V}$$
Generation rate of TEV:
$$g_{T E V}=\frac{V_{T E V, \max }}{1+\left(\frac{\beta_{T F} C_{T F}}{K_{h}}\right)^{n}}$$
Degradation rate of TEV:
$$d_{TEV}=\alpha_{TEV}C_{TEV}$$
So, we get the concentration change of TEV:
$$\frac{d C_{T E V}}{d t}=\frac{V_{T E V, \max }}{1+\left(\frac{\beta_{T F} C_{T F}}{K_{h}}\right)^{n}}-\alpha_{T E V} C_{T E V}$$
For TF:
$$\frac{dC_{TF}}{dt}=g_{TF}-d_{TF}$$
Generation rate of TF:
$$g_{TF}=Const$$
Consider the degradation rate of TF under two conditions. When the concentration of TF is large enough (large than TEV in this model), the degradation rate can be described by M equation. Otherwise, if the concentration is relatively small (smaller than TEV in this model), the degradation rate can be considered to be in direct proportion to the concentration of TF. We get:
$$d_{T F}=\beta_{T E V} \frac{V_{T F, \max } C_{T F}}{K_{m}+C_{T F}}+\alpha_{T F} C_{T F} \quad\left(C_{T F} \geq \beta_{T E V} C_{T E V}\right)$$
$$d_{T F}=\left(k_{d, T F}+\alpha_{T F}\right) V_{T F} \quad\left(C_{T F}<\beta_{T E V} C_{T E V}\right)$$
We get the concentration of TF:
$$\frac{d C_{T F}}{d t}=\operatorname{Const}-\left(\beta_{T E V} \frac{V_{T F, \max } C_{T F}}{K_{m}+C_{T F}}+\alpha_{T F} C_{T F}\right) \quad\left(C_{T F} \geq \beta_{T E V} C_{T E V}\right)$$
$$\frac{d C_{T F}}{d t}=C \text { onst }-\left(k_{d, T F}+\alpha_{T F}\right) C_{T F} \quad\left(C_{T F}<\beta_{T E V} C_{T E V}\right)$$
According to following assumptions, we determine the value of parameters:
1) The generation rate of TEV and TF, and degradation rate that TEV degrades TF are approximately equal. Pick up \( V_{TEV,\ max}=Const=V_{TF,\ max}=1 \).
2) TEV and TF are relatively stable, the self-degradation rates are small. Pick up \( \alpha_{TEV}=\alpha_{TF}=0.05\)
3) Pick up \(K_{h}=K_{m}=1,n=1.5\). Their values make slight influence to the tendency of concentration change of TEV and TF.
4) The transition range of thermosensitive TEV and TF are 40-42℃, the change of efficiency of TEV and TF are linear.
Consider the process that temperature changes from low to high and to low, we simulate the change of concentration TF under conditions of thermosensitive TF and non-thermosensitive TF.
Figure 2. The change of TF’s concentration versus time
Compared with using non-thermosensitive TF, our design has lower TF concentration when temperature is low. Circuit with thermosensitive is a better choice.
2.mf-lon-GOI Subsystem
Parameter used in this model：
| $$\alpha_1$$ | Generation rate of x |
| $$d_{x}$$ | Self-degradation rate of GOI |
| $$\alpha_2$$ | Concentration decrease rate caused by split |
| $$\alpha=\alpha_1+\alpha_2$$ | Sum of above two degradation rates |
| $$V_{max}$$ | The max degradation rate by mflon |
| $$K_m$$ | Coefficient for Michaelis-Menten equation |
| $$P$$ | Concentration of protein |
We introduce mflon to our system to accelerate the degradation of output proteins, consequently, we model the mflon-GOI sub-system to verify and evaluate the acceleration.
The Generation rate of GOI can be seen as zero when the switch is off. Only take the degradation of GOI into consideration.
When there is no mflon, the degradation of GOI depends on the self-degradation and split of bacteria. The concentration change of GOI can be written as:
$$\frac{dP}{dt}=-\alpha P$$
When mflon is imported to the system, assume that the degradation caused by mflon follow the Michaelis-Menten equation, the degradation can be presented as:
$$\frac{d P}{d t}=-\alpha t-\frac{V_{\max } P}{K_{m}+P}$$
According to following assumptions, we determine the value of parameters:
1) The split rate of bacterium is approximate to 30 mins per generation.\(\alpha_2=\frac{ln2}{30} \)
2) The self-degradation rate of GOI is relatively small, can be ignored.\( \alpha\approx 0\)
3) The degradation that mflon degrades GOI is twice as self-degradation.
We simulate the degradation process with system containing mflon gene and not containing mflon gene.
Figure 3. The concentrate of GOI versus time
In system containing mflon, the concentration decreases faster. The time at which concentration decreases to 0.5 in system containing mflon is around 1/3 of the time in system not containing mflon.
3.Determine of Some Parameters
Parameter used in this model：
| $$g_{x}$$ | Generation rate of x |
| $$d_{x}$$ | degradation rate of x |
| $$C_{x}$$ | Concentration rate of x |
| $$\alpha_{x}$$ | Self-degradation rate of x |
| $$\beta_{x}$$ | Relative efficiency of x |
| $$n_i$$ | Hill coefficient |
| $$K_h$$ | K for hill equation |
| $$K_m$$ | K for Michaelis-Menten equation |
| $$K_h$$ | K for hill equation |
| $$V_{g,x,max}$$ | Max generation rate of x |
| $$V_{d,x,max}$$ | Max degradation rate of x |
We design experiments to determine parameters in model, which can help to the describe the gene circuit more precisely.
We first conduct the experiment to determine the response of promoter \(P_{R}\)to TF.
Figure 4. Circuit design of experiment
Under different concentrations of IPTG, and different temperature we measure the fluorescence intensity of GOI.
Consider the concentrate of TF, we get:
$$ \frac{d C_{T F}}{d t}=\frac{V_{g, T F, \max } C_{I P T G}^{n_{1}}}{K_{h 1}+C_{I P T G}^{n_{1}}}+b-\alpha_{T F} C_{T F}$$
When the system reaches a steady state, the concentration of TF changes slightly, the change rate can be seen as 0, i.e.:
$$ \frac{V_{g, T F, \max } C_{I P T G}^{n_{1}}}{K_{h 1}+C_{I P T G}^{n_{1}}}+b-\alpha_{T F} C_{T F}=0 $$
The concentration of TF can be solved:
$$ C_{T F}=\frac{1}{\alpha_{T F}}\left(\frac{V_{g, T F, m a x} C_{I P T G}^{n}}{K_{h 1}+C_{I P T G}^{n T G}}+b\right) $$
Similarly, the concentration of GOI at steady state can be solved:
$$ C_{G O I}=\frac{1}{\alpha_{G O I}}\left(\frac{V_{g, G O I, \max }}{K_{h 2}+C_{T F}^{n_{2}}}\right) $$
Consequently, the concentration of GOI can be written as:
$$ C_{G O I}=\frac{A}{1+B\left(\frac{C_{I P T G}^{n_{1}}}{K_{h 1}+C_{I P T G}^{n 1}}+b\right)^{n_{2}}} $$
Where \( A=\frac{V_{g, G O I, \max }}{\alpha_{G O I} K_{h 2}}\),\(B=\frac{1}{K_{h 2}}\left(\frac{V_{g, T F, m a x}}{\alpha_{T F}}\right)^{n_{2}} \)
We fix \( n_2\) as 1.5, fix \(n_1\) as 1 for the simplification of model.
Choose \(A\) as the max fluorescence intensity from experiment—40000.
Fitting \(K_{h1}\),\(B\) ,and \(b\) in different temperature. Assume that \(B\) changes with temperature because TF is thermosensitive, while \(b\) and \(K_{h1}\) keep the same in different temperature.
Fit the curve with experiment data, and get the parameters as follows:
Table 1 Parameters obtained by fitting
| $$A$$ | $$b$$ | $$K_{h1}$$ |
| 40000 | 0.7 | 0.02623 |
| $$Temperature/℃$$ | $$30$$ | $$34$$ | $$37$$ | $$40$$ |
| $$B$$ | 1.182 | 0.9437 | 4.219 | 3.732 |
The fitting curve in 37℃ is shown:
Figure 5. Fitting curve and experiment data in 37℃
We then conduct another experiment to determine the rate that TEV degrades TF. The circuit designed is shown:
Figure 6. Circuit design of experiment
Under different concentrations of IPTG, and different temperature we measure the fluorescence intensity of GOI.
Consider the concentration of TF:
$$ \frac{d C_{T F}}{d t}=g_{T F}-d_{T F} $$
Assume that the self-degradation rate of TF can be ignored. When the concentration of TF is small, \(k_{d,TF}C_{TF}\) , when the system reaches stable state, the concentration of TF can be solved:
$$C_{TF}=\frac{g_{TF}}{k_{d,TF}}$$
When the concentration of TF is relatively large, \(d_{TF}=\frac{V_{d,TF,max}C_{TF}}{K_m+C_{TF}}\) , when the system reaches stable state, the concentration of TF can be solved:
$$ C_{T F}=K_{m}\left(-1+\frac{V_{d, T F_{,}, m a x}}{V_{d, T F, m a x}-g_{T F}}\right) $$
The concentration of TF can be obtained by retrodict from the fluorescence intensity, and the generation rate can be obtained from last experiment.
The \(k_{d,TF}\) is obtained by \frac{g_{TF}}{C_{TF}}when \(C_{TF}\) , the result is 3.48. When \(g_{TF}\) is large enough, the \(C_{TF}\) tends to be \(-K_m\) . \(K_m\) takes the stable value of \(C_{TF}\) ,which is different in different temperature. \(V_{d,TF,max}\) is obtained by fitting, and the result is -0.5203. The parameters are listed below:
Table 2 Parameters obtained by fitting
| $$V_{d,TF,max}$$ | $$k_{d,TF}$$ |
| -0.523 | 3.5 |
| $$Temperature/℃$$ | $$30$$ | $$34$$ | $$37$$ | $$40$$ |
| $$K$$ | -0.3566 | -1.1118 | -6.4279 | -10.3497 |
Part 2 Drug Distribution
In this part we consider the drug distribution, show that compared to the traditional method, delivery by bacteria can obtain more stable drug concentration in blood.
Drug distribution is relative to the distribution of bacterium and the process of drug delivery. Therefore, we first establish a model for the distribution of bacterium and then consider the delivery of drug.
1.Distribution of Bacterium
In this part we simulate the concentration of bacteria and show that the distribution of bacteria tends to be even after time of several hours.
Parameter used in this model：
| $$n_{i}$$ | Number density of i (scalar field) |
| $$t$$ | Time |
| $$J_{i}$$ | Flow filed of i (vector field) |
| $$\alpha_{i}$$ | Transition rate of i (scalar field) |
| $$D$$ | Diffusivity of bacterium (Const) |
| $$u$$ | Flow velocity of fluid (vector field) |
| $$K_h$$ | K for hill equation |
| $$k$$ | Growth rate of bacterium |
| $$g_i$$ | Generating function of i |
Consider the concentrate change of bacterium after they are released to in the intestinal tract. We use f to represent free bacteria, c to present colonized bacterium.
Assume that the distribution change of bacterium depends only on current bacterium distribution and environment conditions. We get basic equations:
$$\frac{\partial n_{f}}{\partial t}+\nabla \cdot J_{n_{f}}=g_{f}\left(n_{f}, n_{c}\right)$$
$$\frac{\partial n_{c}}{\partial t}=g_{c}\left(n_{f}, n_{c}\right)$$
Consider the diffuse of bacterium and the fluxion of liquid, we get:
$$J_{n_{f}}=-D \nabla n_{f}+n_{f} u$$
Consider the transition and grow of bacterium, we get:
$$g_{f}\left(n_{f}, n_{c}\right)=k n_{f}-\alpha_{f} n_{f}+\alpha_{c} n_{c}$$
$$g_{c}\left(n_{f}, n_{c}\right)=k n_{c}+\alpha_{f} n_{f}-\alpha_{c} n_{c}$$
Substitute to the first two equations, we have:
$$\frac{\partial n_{f}}{\partial t}+u \cdot \nabla n_{f}-D \Delta n_{f}=k n_{f}-\alpha_{f} n_{f}+\alpha_{c} n_{c}$$
$$\frac{\partial n_{c}}{\partial t}=k n_{c}+\alpha_{f} n_{f}-\alpha_{c} n_{c}$$
Here we use that incompressible fluid satisfies the property \( \nabla \cdot u=0 \)
.Now consider the environment of intestinal tract, the intestinal can be treated as one dimension. Then the equation become:
$$\frac{\partial n_{f}}{\partial t}+u \frac{\partial n_{f}}{\partial x}-D \frac{\partial^{2} n_{f}}{\partial x^{2}}=k n_{f}-\alpha_{f} n_{f}+\alpha_{c} n_{c}$$
$$\frac{\partial n_{c}}{\partial t}=k n_{c}+\alpha_{f} n_{f}-\alpha_{c} n_{c}$$
According to following assumptions, we determine the value of parameters:
1) The original concentration of free bacterium is 1.\(n_f(t=0,x=0)=1\).
2) The process that colonized bacterium become free is difficult. \(\alpha_c\approx 0\).
3) The rate that bacterium trans from free to colonized is a const. And we let it be 0.05. \(\alpha_f=0.05\)
4) The diffusivity of bacterium is about \( 3\times10^{-5}cm^{2}/s\). \(D=3\times10^{-5}cm^{2}/s=0.0018xm^2/min\).
5) The vertical of fluid carrying bacterium is about \( 5cm/min \).\(u=5cm/min\).
6) The growth of bacterium in intestinal satisfies logistic curve, and the colonized bacterium growth faster than free bacterium.
$$k_{c}=\frac{\ln (2)}{30} \times \frac{n_{c}\left(K_{c}-n_{c}\right)}{K_{c}}, K_{c}=2$$
$$k_{f}=\frac{\ln (2)}{60} \times \frac{n_{f}\left(K_{f}-n_{f}\right)}{K_{f}}, \quad K_{f}=0.5$$
We simulate the distribution of colonized in 10 hours:
Figure 7. The distribution of bacteria in 600 mins
This figure shows that after around 6-7 hours, the distribution of bacteria tends to be even. Because that the area that capsule can heat is limited to under 10cm, bacteria outside the scope will die of toxin.
2.Drug Delivery
The system is designed to deliver drugs. One of its possible usages is delivering levodopa in the treatment for Parkinson’s disease. Therefore we use this model to see whether the drug delivery of our system can reduce the fluctuation of the drug concentration in plasma, thus will result in a more stable and controllable clinical response.
Before we look into the model, some important notes are listed below:
For the drug in the plasma, as the characteristic time of blood circulation is much shorter than the characteristic time in the pharmacodynamic progress we concern, a one-compartment model will be involved, in which we regard the distribution of the drug in the plasma as homogeneous.
For further consideration, we use the “effect compartment” model to estimate the effect of the drug, which perfectly fits the statistics from the experiment, according to M. Contin(1993) [1].
All the progress involved are linear and can be described by first-order differential equations, by other words, The rate of drug delivery is in direct proportion of the drug amount. This is a fact supported by numerous experiments about levodopa (for example M. Contin(1993) [1]).
The parameters involved are listed below:
| $$f(t)$$ | The drug release rate in the intestine |
| $$C_{\text{max}}$$ | Concentration of drugs in x |
| $$k_{x,y}$$ | Transporting rate of the drug from x to y |
| $$E_{0}$$ | The baseline for effect |
| $$E$$ | The total effect |
| $$E_{max}$$ | The maximum change in effect from baseline |
| $$EC_{50}$$ | the drug concentration at which the change in effect is half maximum |
| $$N$$ | The Hill number |
| $$V_{0}$$ | Bacterium’s generation rate of drugs |
The whole progress consists of three steps. First, the drug is released and gathers in the intestine; second, the drug is absorbed thus enters the plasma; for the last step, the drug enters the “effect compartment”. The following figure (Figure 8) explains the flow of drug. Here we emphasize that the “effect compartment” is only a model and the drug concentration in it cannot be measured directly (thus we can make the maximum of the drug concentration in plasma also the maximum of the “effect compartment”).
Figure 8. the Progress of Drug Delivery of Levodopa
The progress can be described by the following equations.
$$\frac{d C_{\text {intestine }}}{d t}=f(t)-k_{\text {intestine, plasma }} C_{\text {intestine }}$$
$$\frac{d C_{\text {plasma }}}{d t}=k_{\text {intestine, plasma }} C_{\text {intestine }}-k_{\text {plasma }, \varnothing} C_{\text {plasma }}$$
$$\frac{d C_{\text {effect }}}{d t}=k_{\text {plasma, effect }} C_{\text {plasma }}-k_{\text {effect }, \varnothing} C_{\text {effect }}$$
As the system is linear and has the homogeneity of time, we can take \(f(t)=\delta(t)\) ,where the \(\delta(t)\) is the Dirac Delta Function, with which we can get the “basic solution”. Here we use \(C_{\mathrm{x}}^{*}\) to distinguish the “basic solution” from other solutions.
$$C_{\text {intestine }}^{*}(t)=H(t) e^{-k_{\text {intestine,plasma }} t}$$
$$C_{\mathrm{plasma}}^{*}(t)=\frac{k_{\mathrm{intestine}, \mathrm{plasma}}}{k_{\mathrm{intestine}, \mathrm{plasma}}-k_{\mathrm{plasma}, \varnothing}} H(t)\left(-e^{-k_{\mathrm{intestine,plasma}}t}+e^{-k_{\mathrm{plasma}, \varnothing}t}\right)$$
$$C_{\mathrm{effect}}^{*}(t)=\frac{k_{\text {intestine,plasma }} k_{\text {plasma,effect }}}{k_{\text {intestine,plasma}}-k_{\text {plasma,}\varnothing}} H(t)\left(\frac{-e^{-k_{\text {plasma},\varnothing }t}+e^{-k_{\text {effect},\varnothing} t}}{k_{\text {plasma, }\varnothing}-k_{\text {effect},\varnothing}}-\frac{-e^{-k_{\text {intestine,plasma}}t}+e^{-k_{\text {effect},\varnothing}t}}{k_{\text {intestine,plasma }}-k_{\text {effect},\varnothing}}\right)$$
Where the \(H(t)\) is the Heaviside Step Function.
For a more general \(f(t)\), we can use convolution to get the concentration easily.
$$C_{\mathrm{x}}(t)=\left(f * C_{\mathrm{x}}^{*}\right)(t)=\int C_{\mathrm{x}}^{*}\left(t^{\prime}\right) f\left(t-t^{\prime}\right) d t^{\prime}=\int f\left(t^{\prime}\right) C_{\mathrm{x}}^{*}\left(t-t^{\prime}\right) d t^{\prime}$$
Among them, the \(C_{\mathrm{plasma}}\) is a measurable concentration that ralated to side effects, and the \(C_{\mathrm{effect}}\) determines the clinical effect of the drug, which can be illustrated by Hill equation:
$$E=\frac{E_{\max } C_{\text {effect }}^{N}}{E C_{50}^{N}+C_{\text {effect }}^{N}}+E_{0}$$
Some typical value of the parameters involved are listed below (Given by M. Contin(1993) [1])
(Note: In the list we labeled “s” for “stable” patients' statistics, “f” for “fluctuating” patients' statistics, values of parameters differs for these two kind of patients, leading to significant difference in the response to the levodopa (have been described in Project Description), according to M. Contin(1993) [1])
| $$\ln 2 / k_{\text {intestine}, \text { plama }}(\min )$$ | 10 |
| $$\ln 2 / k_{\mathrm{plasma}, \varnothing}(\min )$$ | 50 |
| $$\ln 2 / k_{\mathrm{effect}, \varnothing}(\min )$$ | 150(s) 30(f) |
| $$C_{\text{max}}(\mu \mathrm{g} / \mathrm{ml})$$ | 1.9 (for a standard dose, 100mg) |
| $$E_{0}$$ | 100 |
| $$E_{max}$$ | 50 |
| $$EC_{50}$$ | 0.25(s) 0.65(f) |
| $$N$$ | 2(s) 20(f) |
Here we considered several cases:
Case 1—Oral Delivery of Drugs
This kind of situation is the most common way patients take levodopa. For a single dose, all the drug can be considered as released at a very short period of time, thus as a good approximation, we have \(f(t)=\alpha \delta(t)\) .
Thus
$$C_{\mathrm{plasma}}(t)=\frac{k_{\mathrm{intestine}, \mathrm{plasma}}}{k_{\mathrm{intestine}, \mathrm{plasma}}-k_{\mathrm{plasma}, \varnothing}} H(t)\left(-e^{-k_{\mathrm{intestine,plasma}}t}+e^{-k_{\mathrm{plasma}, \varnothing}t}\right)$$
$$C_{\mathrm{effect}}(t)=\frac{k_{\text {intestine,plasma }} k_{\text {plasma,effect }}}{k_{\text {intestine,plasma}}-k_{\text {plasma,}\varnothing}} H(t)\left(\frac{-e^{-k_{\text {plasma},\varnothing }t}+e^{-k_{\text {effect},\varnothing} t}}{k_{\text {plasma, }\varnothing}-k_{\text {effect},\varnothing}}-\frac{-e^{-k_{\text {intestine,plasma}}t}+e^{-k_{\text {effect},\varnothing}t}}{k_{\text {intestine,plasma }}-k_{\text {effect},\varnothing}}\right)$$
To meet with the value of \(C_{\text{max}}\) , we have \(\alpha=2 \cdot 5^{\frac{1}{4}} \mu \mathrm{g} / \mathrm{ml} \approx 3.0 \mu \mathrm{g} / \mathrm{ml}\), similarly we can get the value of \( k_{ \text{plasma,effect}} \).
With the parameters’ value given above, we can present the result in the form of graphs.
(Note: The experiemtal data points invloved in Figure 9 are from Urszula Adamiak et al., scales modified.)
Figure 9. Drug Concentration Versus Time, Unit Dose Taken Orally
Figure 10. Clinical Effect Versus Time, Unit Dose or Double Doses Taken Orally
We can see that especially for fluctuate patients, the clinical effect drops quickly back to baseline in about 3 hours after taking the dose. Adding the size of the dose makes no great changes to the situation (see Figure 10), therefore the only solution left to prevent the “wearing-off” phenomena is to take levodopa regularly. The frequent intake of medicine is obviously a great inconvenience for the patient. Moreover, if the patient does not take levodopa frequent enough, or occasionally forget to take it, the clinical effect for fluctuate patients will go up and down like a roller-coaster (see Figure 11).
Figure 11. Clinical Effect for 24h, Unit Dose Taken Orally Every 6h
Nowadays many levodopa medicine are designed to be sustained-released, and we also take it into consideration. In the sustained-release situation, we assume \(f(t)=\frac{1}{\sqrt{2 \pi \sigma^{2}}} e^{-\frac{t^{2}}{2 \sigma^{2}}}\), and take \(\sigma=60 \mathrm{min}\). By numerical calculation we get:
Figure 12. Drug Concentration Versus Time, for Sustained-released Medicine
Figure 13. Clinical Effect Versus Time, the Difference Between None-sustained-released and Sustained-released Medicine
We see that sustained-released medicine is better than none-sustained-released once, however the improvement is not that significant.
Case 2—Bacterial Delivery of Drugs
This case is related to the system we designed. According to the result of last model, distribution of bacterium can be considered as even, so the rate of drug release can be considered as a constant, which means \(f(t)=V_{0}\).
Therefore it is easy to have:
$$C_{\mathrm{plasma}}(t)=\left(f * C_{\mathrm{plasma}}^{*}\right)(t)=\frac{V_{0}}{k_{\mathrm{plasma}, \varnothing}}$$
$$C_{\text {effect }}(t)=\left(f * C_{\text {effect }}^{*}\right)(t)=\frac{k_{\text {plasma }, \text { effect }} V_{0}}{k_{\text {effect }, \varnothing} k_{\text {plasma }, \varnothing}}$$
The relationship between concentration, effect and drug produce rate are illustrated by the graphs below.
(Note: In this situation, the concentration of the drug as well as the clinical effect does not change with time, therefore the figures below are not about the changes of certain variables versus time.)
Figure 14. Drug Concentration Versus Drug Produce Rate, Stable Delivery of Medicine by Bacterium
Figure 15. Clinical Effect Versus Drug Produce Rate, Stable Delivery of Medicine by Bacterium
Therefore by adjusting the amount of bacteria and the productivity of bacteria (like if we adjust \(V_{0} \geq 0.0025 \mu \mathrm{g} /(\mathrm{ml} \cdot \mathrm{min})\), which is much less than the maximum productivity of the bacteria (over \(0.01 \mu \mathrm{g} /(\mathrm{ml} \cdot \mathrm{min})\), see Project Description)), we can achieve an ideal and stable concentration of levodopa, so do the clinical effect.
Now we consider the occasion that the production of drug is influenced by environment, and it fluctuates. We will see that our system is still robust enough to hold a stable effect even for the fluctuate patients. We take \(f(t)=\frac{(1-\cos \omega t)}{2} V_{0}\) (it is too fluctuate to be a real one, thus can show more about the robustness of our system), and simulates by numerical calculation. Here are the results.
Figure 16. Drug Concentration Versus Time, Unstable Delivery of Medicine by Bacterium
Figure 17. Clinical Effect Versus Time, Unstable Delivery of Medicine by Bacterium
When \(V_{0}=0.003\mu \mathrm{g} /(\mathrm{ml} \cdot \mathrm{min})\), the clinical effect for the stable patients has already reached the maximum, but for the fluctuate patients it is very low (at about baseline value). While \(V_{0}=0.005\mu \mathrm{g} /(\mathrm{ml} \cdot \mathrm{min})\), the curve of the clinical effect is becoming stable for fluctuate patients. System’s high ability of staying stable is in fact a consequence of the convolution effect, making our system more advanced than the oral delivery of levodopa.
a More Detailed Comparison of the Way of Levodopa Delivery
We may first find out the origin of the drawbacks of oral delivery of drugs. We again look at the Hill equation that describes the clinical effect
$$E=\frac{E_{\max } C_{\text {effect }}^{N}}{E C_{50}^{N}+C_{\text {effect }}^{N}}+E_{0}$$
By its dimension and definition, we can define \(EC_{50}\) to be the “efficient concentration”— a threshold value below which there is little clinical effect while above which the clinical effect will not increase much. This especially fit for high Hill number situation, like for fluctuate patients (\(N=20\)).
For oral delivery, we can see that the peek of the drug concentration is far higher than the efficient concentration (see Figure 18), which can bring no significant increase of clinical effect. Moreover, for fluctuate patients, the concentration drops too fast so that the efficient time is short.
Figure 18. Drug Concentration Versus Time, Whith "efficient concentration"
Obviously, the best delivery method should provide a stable drug concentration just above the efficient concentration. The bacterial delivery perfectly meets this need.
Here we made a comparison between bacterial delivery and oral delivery, by comparing the clinical effect versus time, and by comparing the drug concentration in plasma, which is related to possible side effects.
(Note: We considered a time period of 24h. For oral delivery, the patient take a unit dose of levodopa every 6h. For bacterial delivery, we considered both stable and unstable situation, for both situation, the mean drug produce rate is \(V_{0}=0.003\mu \mathrm{g} /(\mathrm{ml} \cdot \mathrm{min})\) for stable patients, and \(V_{0}=0.005\mu \mathrm{g} /(\mathrm{ml} \cdot \mathrm{min})\) for fluctuate patients.)
Figure 19. Drug Concentration Versus Time by Different Methods, for Stable Patients
Figure 20. Clinical Effect Versus Time by Different Methods, for Stable Patients (the origin of y axis is 140)
Figure 21. Drug Concentration Versus Time by Different Methods, for Fluctuate Patients
Figure 22. Clinical Effect Versus Time by Different Methods, for Fluctuate Patients
We can see from the figures that the bacterial delivery generally has a high-and-stable clinical effect, while having a much lower drug concentration in plasma, thus will bring much lower side effects. Its advantage over oral delivery is especially significant for fluctuate patients.
Part 3 Hardware
1.Heating
In this part, we simulate the process of heat conducts from our capsule to intestine, and get the temperature distribution of intestine, combined with result of the experiment, to verify the feasibility.
Parameter used in this model：
| $$\rho$$ | Density of tissue |
| $$C$$ | degradation rate of x |
| $$C_{x}$$ | Specific heat capacity of tissue |
| $$K$$ | Thermal conductivity of tissue |
| $$K_{b}$$ | Increase of thermal conductivity |
| $$Q_{b}$$ | Heat exchange caused by hemoperfusion |
| $$\omega_b$$ | Hemoperfusion |
| $$\rho_b$$ | Density of blood |
| $$C_b$$ | Specific heat capacity of blood |
| $$Q_{i}$$ | Other heat |
| $$T$$ | Temperature of tissue |
| $$T_{b}$$ | Temperature of arterial blood |
According the heat conduction applying to biology tissue by Penns[2], we get:
$$ \rho C \frac{\partial T}{\partial t}=\nabla \cdot\left(\left(K+K_{b}\right) \nabla T\right)+\omega_{b} \rho_{b} C_{b}\left(T_{b}-T\right)+Q_{i} $$
This figure by Gautherie(1969)[3] shows the relationship of temperature and the conductivity \(K+K_b\).And Hurlburt and Chato(1987)[4] fitting the curve by the following function.
$$ K+K_{b}=4.82-4.44833 \times 1.00075^{-1.575^{T-25}} $$
Figure 23. The fitting curve by Hurlburt and Gautherie
The intestine can be treated as one dimension. We get:
$$ \rho C \frac{\partial T}{\partial t}=\left(K+K_{b}\right) \frac{\partial^{2}}{\partial x^{2}} T+Q_{b}+Q_{m} $$
$$ \rho C \frac{\partial T}{\partial t}=\left(K+K_{b}\right) \frac{\partial^{2}}{\partial x^{2}} T+\omega_{b} \rho_{b} C_{b}\left(T_{b}-T\right)+Q_{i} $$
$$ K+K_{b}=4.82-4.44833 \times 1.00075^{-1.575^{T-25}} $$
According to following assumptions, we determine the value of parameters:
1) The specific heat capacity of blood and tissue is approximately equal to the heat capacity of water.\(C_{b} \approx C \approx C_{\text {water}}=4200 \mathrm{J} /\left(\mathrm{kg} \cdot \mathrm{C}^{\circ}\right)\).
2) The density of blood and tissue is approximately equal to the density of water. \(\rho_{b} \approx \rho \approx \rho_{\text {water }}=1000 \mathrm{kg} / \mathrm{m}^{3}\).
3) The hemoperfusion is about \( 0.0005 \mathrm{ml} /(\mathrm{s} \cdot \mathrm{ml}) \) .\( \omega_{b} \approx 0.0005 \mathrm{ml} /(\mathrm{s} \cdot \mathrm{ml})\).
4) The temperature of arterial blood is about \(37 C^{\circ}\) .\(T_{b} \approx 37 C^{\circ}\).
The boundary conditions:
$$ T(0, t)=50 C^{\circ} $$
$$ T(0.08, t)=37 C^{\circ} $$
And initial conditions:
$$ T(x, 0)=37 C^{\circ} $$
Using above boundary conditions and initial conditions, we simulate the temperature distribution of intestine in 1000 seconds:
Figure 24. The temperature of heating area in 1000s
The heating range is about 4-5 centimeters in 1000 seconds, which is not too large, effecting colonized bacteria in intestine, not too small, causing little drug effect or target miss, proving that the design of heating capsule is feasible.
2.Location
The heating capsule using external device with magnetic sensors to locate. The location of capsule is a optimization problem, and here is the description of location algorithm[5].
a) Input and output
Input: The 3-dimention coordinates of n sensors.
$$(x_1,y_1,z_1);(x_2,y_2,z_2);\cdots; (x_n,y_n,z_n)$$
The magnetic field intensity at positions of every sensors.
$$(B_{1_x},B_{1_y},B_{1_z});(B_{2_x},B_{2_y},B_{2_z});\cdots(B_{n_x},B_{n_y},B_{n_z})$$
Output:
The 3-dimention coordinates of magnet.
b) Modeling method--optimization method
Assume that the magnet can be seen as a magnetic dipole, then suppose a position of a magnetic sensor to be \((x_l,y_l,z_l)\) , the position of magnet to be\((x_d,y_d,z_d)\), and the posture of the magnet can be described by the cosine of angles between the line connecting N and S poles and axis--\((m,n,p)\) .
Then the magnetic field density is supposed to be \((B_x,B_y,B_z)\) where:
$$ B_{x}=B_{T}\left\{\frac{3\left(x_{l}-x_{d}\right)\left\{m\left(x_{l}-x_{d}\right)+n\left(y_{l}-y_{d}\right)+p\left(z_{l}-z_{d}\right)\right\}}{R_{l}^{5}}-\frac{m}{R_{l}^{3}}\right\} $$
$$ B_{y}=B_{T}\left\{\frac{3\left(y_{l}-y_{d}\right)\left\{m\left(x_{l}-x_{d}\right)+n\left(y_{l}-y_{d}\right)+p\left(z_{l}-z_{d}\right)\right\}}{R_{l}^{5}}-\frac{n}{R_{l}^{3}}\right\} $$
$$ B_{z}=B_{T}\left\{\frac{3\left(z_{l}-z_{d}\right)\left\{m\left(x_{l}-x_{d}\right)+n\left(y_{l}-y_{d}\right)+p\left(z_{l}-z_{d}\right)\right\}}{R_{l}^{5}}-\frac{p}{R_{l}^{3}}\right\} $$
$$ R_{l}=\sqrt{\left(x_{l}-x_{d}\right)^{2}+\left(y_{l}-y_{d}\right)^{2}+\left(z_{l}-z_{d}\right)^{2}} $$
$$B_T=\frac{\mu _r \mu_0M_T}{4\pi}$$
\(M_T\) represents the Magnetization.
Supposed the magnetic field density measured by magnet sensor to be \((B_{lx},B_{ly},B_{lz})\)
The error function for one sensor is defined as:
$$E=(B_x-B_{lx})^2+(B_y-B_{ly})^2+(B_z-B_{lz})^2$$
Given the coordinates of n sensors and the magnetic field density at them, the error function can be written as:
$$ E=\sum_{i=1}^{n} E_{i}=\sum_{i=1}^{n}\left(\left(B_{x_{i}}-B_{l x_{i}}\right)^{2}+\left(B_{y_{i}}-B_{l y_{j}}\right)^{2}+\left(B_{z_{i}}-B_{l z_{i}}\right)^{2}\right) $$
The location algorithm is to find a set of parameters to minimize E.The coordinates obtained by this algorithm are considered as the real position of magnet.
c) Calculation method--Levenberg-Marquardt Algorithm
The Levenberg-Marquardt Algorithm is a calculation method for optimization method. Its input is a vector function \(f:R^m\to R^n\) with \(n\ge m\) , a measurement vector for \(x \in R^{n}\) and an initial parameters estimate \(p_{0}\in R_{m}\)
Its output is a vector \(p_{+}\in R^m\) minimizing \(||x-f(p)||^2\)
The algorithm can be described in pseudocode as follows:
\(k:=0;v:=2;p:=p_0;\)
\(A:=J^{T}J;\epsilon:=x-f(p);g:=J^T\epsilon_p;\)
stop\(:=(||g||_{\infty}\le \epsilon_1);\mu:=\tau*max_{i=1,2,\dots,m}(A_{ii});\)
while(not stop)and\((k\le k_{max})\)
\(k:=k+1;\)
repeat
solve\((A+\mu I ) \delta _p=g\)
if(\( ||\delta_p||\le\epsilon_2||p||\))
stop\(:=\)true;
else
\(p_{new}:=p+\delta_{p};\)
\(\rho:=(||\epsilon_p||^2-||x-f(p_{new})||^2)/(\delta_p^T(\mu\delta_p+g))\)
if(\(\rho>0\))
\(p=p_{new};\)
\(A:=J^{T}J;\epsilon_{p}:=x-f(p);g:=J^{T}\epsilon_{p};\)
stop:=\((||g||_{\infty}\le \epsilon_1)\)or\((||\epsilon_p||^2\le\epsilon_3)\)
\(mu:=mu*max(\frac{1}{3},1-(2\rho-1)^3);v:=2\)
else\(\mu:=\mu*v;v:=2\times v\)
endif
endif
until\((\rho>0)or(stop)\)
endwhile
References
[1]M. Contin, R. Riva, P. Martinelli, et al.(1993), Pharmacodynamic modeling of oral levodopa: Clinical application in Parkinson's disease, Neurology 1993;43;367.
[2]H.H.Pennes (1948), Analysis of Tissue and Arterial Temperatures in the Resting Human Forearm, Journal of Applied Physiology. Vol.1, pp.93-122.
[3]Gautherie M (1969) Etude par thermometrie infrarouge des properties thermiques de tissues humains "in vivo". Influence de la temperature et de la vascularization. Rev Fr Etud Clin Bioi 14:885-901
[4] Hurlburt E, Chato JC (1987) Computer aided modeling of tissue temperatures using finite difference methods. In: BT Chao symposium volume, Department of Mechanical and Industrial Engineering, University of Illinois, Urbana
[5]K.Madsen, H.B.Nielsen, O.Tingleff (2004), Methods for non-linear least squares problems.
