Towards a computational prediction for the tumor selective accumulation of paramagnetic nanoparticles in retinoblastoma cells

Johansen Rj; Bukhvostov Aa; Ermakov Kv; Kuznetsov Da

TOWARDS A COMPUTATIONAL PREDICTION FOR THE TUMOR SELECTIVE ACCUMULATION OF PARAMAGNETIC NANOPARTICLES IN RETINOBLASTOMA CELLS

Johansen RJ1>2, Bukhvostov AA3, Ermakov KV3, *Kuznetsov DA2,3 23

1 Department of Mathematics and Computer Science, University of Southern Denmark, Odense, DK-5230, Denmark

2 N. N. Semenov Institute of Chemical Physics, Russian Academy of Sciences, Moscow

3 Department of Medical Nanobiotechnologies, N. I. Pirogov Russian National Research Medical University, Moscow

Retinoblastoma is a malignant growth affecting retina. An original combination of modified Non-Markov and Gompertzian computational approaches is proven of being a reliable tool for prediction of tumor selective accumulation of the bivalent metal isotopes (25Mg, 43Ca, 60Co, 67Zn, ...) — releasing nanoparticles in human retinoblastoma cells. This mathematical model operates with a starting point of the discriminative drug uptake caused by a gap-like distinction between the neighboring malignant and normal cell proliferation rates. This takes into account both pharmacokinetic and pharmacodynamic peculiarities of PMC16, fullerene-C60 based nanoparticles, known for their unique capabilities for a cancer-targeted delivery of paramagnetic metal isotopes followed by an essential chemotherapeutic effect. Being dependent on a tumor growth rate but not on the neoplasm steady state mass, a randomized level of drug accumulation in retinoblastoma cells has been formalized as a predictive paradigm suitable to optimize an ongoing PMC16 preclinical research.

Keywords: retinoblastoma, paramagnetic cytostatics, nanocationites, tumor selective nanoparticles uptake, drug accumulation mathematical model

Acknowledgments: this work was performed due to an exceptional technical assistance kindly provided by Erasmus-Plus DK06811/2020 Program associates affiliated with the Southern Denmark University at Odense, Denmark, and, most specifically, by Ms. Patricia Wladycziewski, SDU Erasmus chief supervising officer.

[23 Correspondence should be addressed: Dmitry A. Kuznetsov Ostrovityanova 1, Moscow, 117997; [email protected]

Received: 27.06.2018 Accepted: 18.08.2018

DOI: 10.24075/brsmu.2018.078

МАТЕМАТИЧЕСКОЕ ПРОГНОЗИРОВАНИЕ ПАРАМЕТРОВ ОПУХОЛЬ-СЕЛЕКТИВНОГО НАКОПЛЕНИЯ ПАРАМАГНИТНЫХ НАНОЧАСТИЦ КЛЕТКАМИ РЕТИНОБЛАСТОМЫ

Р. Дж. Йохансен1'2, А. А. Бухвостов3, К. В. Ермаков3, Д. А. Кузнецов2^

1 Кафедра математики и компьютерных наук, Университет Южной Дании, Оденсе, Дания

2 Институт химической физики имени Н. Н. Семенова РАН, Москва

3 Кафедра медицинских нанобиотехнологий, Российский национальный исследовательский медицинский университет имени Н. И. Пирогова, Москва

Ретинобластома — злокачественное новообразование, поражающее сетчатку глаза. Целью работы было разработать вычислительный подход к прогнозированию опухоль-специфического накопления наночастиц, высвобождающих катионы изотопов двухвалентных металлов (25Mg, 43Ca, 60Co, 67Zn, ...) в клетках ретинобластомы человека. Предложена математическая модель, основанная на применении уравнения Гомперца и оригинальной версии немарковской популяционной динамики. Она основана на факте ярко выраженного дискриминационного распределения препарата между злокачественными и «соседствующими» с ними нормальными клетками и различиях в параметрах их клеточных циклов. Учтены как фармакокинетические, так и фармакодинамические особенности наночастиц PMC16 — порфирин-производных фуллерена C60, известных благодаря их уникальным возможностям в отношении направленной доставки парамагнитных изотопов металлов в раковые клетки, сопровождающейся существенным химиотерапевтическим эффектом. Демонстрируя зависимость от скорости роста опухоли, но не от ее массы в стационарной фазе, рандомизированный уровень накопления препарата в клетках ретинобластомы формализован как ценный в прогностическом отношении расчетный метод, пригодный для оптимизации проводимых в настоящее время доклинических исследований катионообменных наночастиц PMC16.

Ключевые слова: ретинобластома, парамагнитные цитостатики, нанокатиониты, опухоль-селективное накопление наночастиц, математическая модель накопления препарата

Благодарности: эта работа была бы невозможна без помощи в технических вопросах, которую нам оказали сотрудники Университета Южной Дании в Оденсе в рамках программы Erasmus-Plus DK06811/2020. Авторы выражают особую благодарность Патриции Владычевской, главному руководителю SDS Erasmus.

[X] Для корреспонденции: Дмитрий Анатольевич Кузнецов ул. Островитянова, д. 1, г Москва, 117997; [email protected]

Статья получена: 27.06.2018 Статья принята к печати: 18.08.2018

DOI: 10.24075/vrgmu.2018.078

Human retinoblastoma (RB) is found to be very sensitive to some metal paramagnetic isotopes due their ability to promote a so called magnetic isotope effects which, in turn, promotes a sharp inhibition of DNA repair in malignant cells along with a formation of shorted, and hence DNA repair inconsistent, DNA sequences [1-4]. This might be taken as a "hopeful pullout" for coming up with a new element in RB chemotherapy based on administration of 25Mg2+, 43Ca2+, 60Co2+, 67Zn2+ carrying/ releasing nanoparticles (NPs) once the RB cell does indeed overexpresses the DNA Polymerase Beta, a target enzyme for the nuclear spin selective DNA repair [1, 5, 6].

These complexes of paramagnetic isotopes with PMC16 (Fig. 1), a peculiar type of amphiphilic low-toxic NPs, were in fact deliberately developed to face a requirement for ion transporter applicable in both tumor cell targeting and a subsequent intracellular controlled drug release [1, 5]. As a sign of such paramagnetic impacts, a significant decrease of proliferation rates has been observed in Y79 and WERI-RB-1 retinoblastoma cell strains [2-4].

According to PubMed statistics, the amount of publications on nanoparticles (NPs) for a passive targeted drug delivery has been increased in the past 15 years from about 40 (year 2000) up to nearly 1,800 (2015) taking the solid tumors research only [7]. As per the PMC16 passive targeting which would presumably take place in RB engaging preclinical studies, a tumor selectivity of anticipated NPs uptake looks rather obscure and unpredictable owing to a number of the RB-marking epigenetic factors [8-10].

A reliable prediction on the rate and extent of NP (PMC16) — RB selective accumulation would be no doubt a sort of beneficial supplement to anti-RB chemotherapeutic strategies proposed for a preclinical trial program. This work is an attempt to solve this task by employing a certain arsenal of mathematical modeling tools.

Noteworthy, an autonomous trend of computational approach has already made an essential contribution to preclinical and clinical trial scenaria in oncology and related areas [11-17].

METHODS

To proceed the simulation data, the most common drug (NP) — cell distribution and the RB/RT cell proliferation patterns (Table) [8, 9, 15, 18] were treated using a Sigma QXL600 software algorithm in HP9107 (Hewlett-Packard, Inc; USA) and Olivetti Riccetta SL110 (Ing. C. Olivetti & Co.; Италия) analytical units adopting a slightly modified Penman-Dalbreaux probabilistic approximation technique [14, 16] to harmonize the output with the population dynamics platforms based on both non-Markov [12, 14] and Gompertz [11] equation systems.

RESULTS

Drug toxicity to normal tissues and the emergence of drug-resistance along with a tumor selectivity in drug (NP) targeting/ accumulation processes are no doubt the major limiting barriers on a path to chemotherapy of cancer [5, 6, 8]. A computational modeling of cell population dynamics in harshly varying object-surrounding environment could be applied to chemotherapeutic paradigm [7, 14, 17, 18]. In several cases, this approach might make a difference for improving responsiveness to the phase-specific drugs (NPs) taking into account their non-discriminative, vector-free ("passive"), distribution within a cell pool consisting of neighboring slow and fast proliferating populations.

Non-Markov population dynamics

The dynamics under various drug regimens of populations that differ in life-cycle parameters is simulated using a computer model whose simplest form is given in:

x(t) = X x (t- t)[1 - D(t)],

where x(t) is population density at time t, A is the cell birth rate, t is generation time, and D(t) denotes the environmental process, so that D(t) > 0 corresponds to the occurrence of effective concentration of the drug in the system. Using this model the elimination time of malignant population (7m) and that of the limiting host population (Th) were estimated, and the elimination coefficient, Z, measuring the treatment efficacy, was calculated according to:

Z = 1 -

T IT..

m h

The treatment efficacy is a nonmonotonic function of the relation between the cell generation time and the period of drug administration, with maximal occurring when the limiting host cell cycle length is a multiple of the chemotherapeutic period. Analytical results further show that in fully periodic systems elimination time, T, is given for t > 5 > t/2:

T = tw/|t - (5 + w)| .

Here, 5 is the duration of the period in which the drug effective, and w is the period in which the drug dosage is below efficiency. The point, t = 5 + w, is a singular point with T being infinite.

This makes possible to assume that a classical non-Markovian model of population dynamics [12, 14] is indeed an appropriate tool to simulate the NP (PMC16) distribution between malignant (fast expanding compartment) and the hostile normal cell (slow expanding compartment) pools. The above mentioned amphiphilic pharmacophore (PMC16, cyclohexyl(C60)porphyne-based bivalent metal isotopes nanocarrier; Fig. 1) is a suitable probe for our non-Markovian simulation since this type of NPs was found capable to manifest a clear and sharp cytostatic mode in acute myeloblast leukemia and RB cell cultures [1, 5, 6].

A two-compartment model we proposed is fitted to the following non-Markovian compatible pharmacokinetics data with both inter-specimen and randomized effects on CL, V, Q, and V2 corrected to an error best described the pattern of residual error [12, 13, 16]. So this our model works out for both PMC16 tumor uptake selectivity (fast proliferation caused

Fig. 1. Structure of PMC16 (cyclohexyl(C60)porphyrin), Me2+ — carrying and releasing nanoparticles with the marked membranotropic/amphiphilic properties [1]

Table 1. Population turnover in Y79 and WERI-RB-1 cell Unes

Parameter Meaning T, hr Ref

TG1 Duration of G phase 8.0 [9, 15]

TS Duration of S phase 7.5 [8, 15]

Tg2m Duration of G2M phase 2.0 [8, 18]

Tg0 Duration of G0 phase 3.0 [9, 15]

TApoptosis Duration of the apoptotic phase 4.0 [8, 9]

phenomenon) and a routine pharmacokinetic key points prediction.

Inter-specimen covariant models

A. Pharmacokinetic model (a - k2l)

C = D/V •

exp - (a • t) + ^ ß) • exp - (ß • t)

(a - ß) r v (a - ß)

B. Non-vectoral covariative model

CL. = [03 • OCC1 + e4(WT - 75)] • exp (nCL. )

V. = [ei - (GFR - 80) • e2] • exp (^ )

k12. = (e5) • exp (nti2y )

C. Population dynamics model

, cl , CL , Kr^r'exP- (-rrO

c =-

D

j ei . occ1 + e2 • (wt - 75) • exp (nv . )

CL

ß - -V-1

+ k2l - kl2 CL

-vL- ß

• exp (Ej

• exp - (k12/)

D. Population parameters Gompertzian model

e1 = 19.5

e2= 0.198

a,^

- V • (kvp) =£p .

The computer program starts from an initial distribution of neighboring RB and RT cells in each state {a, 4>}. The compupations are performed using a splitting technique. We run the cell cycle model for one time-step dt, then retrieve new values for n and compute P This drives to a system:

dn

—^ + V • (un ) = 0 dt

dn

dt

■ + u

• vn = ( Z P , ) n

Applied to the cell division cycle key patterns (Table) represented as a non-Markov population dynamics model organized in a merry-go-round of subpopulations biologically identified as phase (G1, S, G2 and M), this might be re-formalized as:

dn(t, x) dn(t, x) „ s „ s

' ■ + —•—— + d (t, x)n (t, x) + K . l(t, x)n.(t, x) = 0

The Gompertz equation based models were already used to describe cancer growth dynamics [7, 11, 14], these formalisms have been also applied to optimize some therapeutic strategies dealing with antiangiogenic [11, 12] and radiation treatment [11, 13].

The model is fully deterministic. Cell cycle phases durations t have been discretized in several elementary age intervals a E {1,..., N+ } where N+ is an integer such as t+ = dt • N . Here dt is the time step of the cell cycle model. The cell density n at age a in phase 4> is governed by:

dn

-^f* + V • (vn ) = P .

In this equation, 4> E {G1, S, G2 M, G0, Apoptosis and a E {1,.,N^}. Pa + is the cell density proliferation term in phase at age retrieved from the cell cycle model. In these simulations, the intracellular and extracellular conditions were identified for cells at the end of G1 phase.

Furthermore, noting that X „ n „ is constant, so we can sum to obtain an expression for the pressure field of the form:

dt dt

nT+1(t, 0) = /X^a x)n.(t, x)dx

ni(t, 0) = 2 J K^(t, x)n(t, x)dx .

along with the initial conditions (n . = 0 )1<.<(.

Cell death rates in phases are noted d^ and transition rates between phases, assumed to be time-periodic. K.^.+1. Phase i (1 < i < () may be one of the classical four G1, S, G2 and M, but also an aggregated phase such as S-G2, or even a single proliferating phase G1-S-G2-M, or, on the contrary, a subdivision inside a phase, e.g., pre- or post-restriction point in G1; the equation describes the evolution of the densities n. (t, x) of cells having age x at time t in phase i.

The above stated two systems that represent two neighboring, fast and slow growing, cell populations are physically apart from each other. Hence, in this system of equations, function g, which represents anti-tumor drug efficacy, is assumed, as is function f for cytotoxicity:

g(D, t) = ll(1 + cos(2n

( t - 9B) 24

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

) )

(DY,B)

DYB + DYB

50

whereas A, v, en, a, B , H, m„, Y„, D,. are positive constants,

' ' D' ' max' ' tb> b 50 " '

identified on tumor growth curves or from literature data [8, 15, 18], or else estimated.

1.0--

^0.5--

■■1.0

■■0.5 H

Fig. 2. The NP uptake selectivity prediction in a complete accessibility of intracellular ligands. P — [NP] uptaken, units per cell; P0 — intracellular initial concentration of NP-ligands; KS — Gompertz equation vectoral K; K — an NP uptake steady state constant; K2 — an efficient constant of saturation of a cellular ligand pool at [NP] ^ 0.5 P

The difference of behaviors between these two populations of cells (RB-RT pair) with respect to drug response is coded as qA - = 13 hours.

Turning back to the roots, a damped harmonic approximation stands for healthy (RT) cell population dynamics:

f_ -ИС + scp

where

^ = {-a - /( C, t)} Z - pA + Y

dA_z z

dt '

/(C t) = F (1 + cos(2n ))

с y,A

rY,A , f Y,A C50 + C

aB ln (jB- - g(D, t)B).

Clearly, this is nothing but the way to represent exchanges with quiescent population in a still linear model which normally means to exclude feedback from quiescence to proliferation, considering quiescence only as a sideway expansion cell tank:

-|-pft + pft x) + In + K(x)}p(t, x) = 0 p(t, x = 0) = 2(1 - f) f K(e) p(t, S) de

S>0

p(t, x = 0) = p0 (x) -dd-Q(t) = 2 f I K(e) p(t, e) de - vQ(t)

dt e>o

Q(0) = Qo

and A, m, £c, a, P, Y, Zeq, F, YA, C50 are positive constants, which, again, were identified on tumor growth curves or from literature data [7, 10, 15], or else estimated.

These equations represent drug diffusion and elimination by first order pharmacokinetics for concentrations in the plasmatic and target cell compartments (P and c), from infusion in the general circulation according to the instantaneous drug delivery flow i (t) ( 0 representing a "tap on-tap off" function), and health tissue homeostasis by a linear system showing a stable focus at Z , A = P-1 (Y - aZ ), perturbed by the drug cytotoxicity function

eq eq eq

which comes to strengthen the natural self-regulation coefficient a.

So our model, as derived from a Gompertz equations row, is completely adequate to the tumor cell population dynamics:

= -\p + -vu(0

j t-\ dist

§= -VD + Sdp B

To emphasize a perspective proclaimed, this our model is to reveal the action of a cytostatic drug enhancing the way out of proliferating cells with density p(t, x) to quiescent cells with density Q (t), the drug target here isfrate of escape at mitosis towards the siding phase Q, f to be enhanced by a cytostatic drug.

DISCUSSION

Tumor selective NP uptake. Probability and prediction

As seen from above, a probabilistic prognosis for the RB-selective NP uptake relates predominantly on a ratio between malignant and normal cell proliferation rates while the mass of a cancer tissue per se (amount of RB cells) might be practically neglected (Fig. 2). This derives from a predictive cell response to a rapidly in situ diffusing probe (PMC16) once these Me2+-nanocarriers arrive to the RB/RT frontier area. In this stochastic scenario, however, a cellular lattice is nothing but a peculiarity reflecting the target cell specific energy landscape [16, 19] which makes the drug trapping probability dependent on the EL motion and, therefore, on expanding dynamics of a most rapidly growing compartment within a given RB/RT pair (Fig. 3).

Z

Fig. 3. Probabilistic model for NPs distribution between RB and the neighboring RT cells as a function of the discriminative cell cycle turnover. Z-ellmination coefficient for malignant and the RB-surrounding normal cells (RT) estimated for the drug efficiency duration time (o) and the drug-free cell functioning time interval (a>). o is normally distributed within a variation range of o = o/10, where remain constant while the inner rate of the "newborn" cell appearance is A = 2, for a starting population size x(0) = 5

A symbolic blue-red shift in Figures 2 and 3 marks a trend to predominant accumulation of NPs in the most faster expanding cell tank, RT.

So turning back to a background proceeding probabilistic approach [12, 14], a tumor selective accumulation of the PMC16-specific probe become predictable due to enormous difference between RB and RT growth rates [8-10, 12]. This allows a rate-discriminative RT-PMC16 uptake described by our model (Fig. 2) working in accordance with:

Ac=Kd[tgaRB/(tgaRB - tgaj],

where Kd is an Arnauld-Pitot disclaimer approximation constant [16,19].

Paul Ehrlich's "magic bullet": dream or nightmare?

Meaning the end of a long lasting post-Virchowian era, a truly prophetic outlook stated by Paul Ehrlich back in 1908-1913, now well-known as a hope for an infamous magic bullet in cellular pathology and pharmacology [8, 17], has been eventually adopted within a contemporary drug targeted delivery concept [5, 7, 19]. The latter requires a broad variety of nanodevices, all sorts of the magic bullets, designed to conduct both the towards-a-target navigated drug transfer and a consequent in situ controlled drug release [1, 5, 17].

However, a new unclear horizon appears straight in front of a marksman equipped with the magic bullet loaded "cartridge". Suppose a reasonable amount of the active drug molecules or ions have reached the tumor location border due to a perfect delivery performed by some nanocarrier.

A good shot with the bullet though. Then a tumor selective intralization of a drug becomes a rather obscure step in a whole pharmacokinetic scenario. Indeed, no matter how precise the bullet's trajectory is, a further distribution of NPs between the neighboring malignant and normal tissue compartments is the case.

That's why a predictive model we proposed is in fact a promising tool the one might need to come up with the trustworthy path/dose/exposition plan to follow and a strategic roadmap to observe upfront of experiment.

After all, a direct Schlemm channel drug influx and/or the intraocular administration paths, often applicable to RB particular case [8], would make this approach not only possible but even preferable as well. A numerous holistic impacts, ineluctable when the routine parenteral administration paths involved, should be therefore minimized or merely neglected as long as the RB chemotherapy is in a focus.

CONCLUSIONS

A mathematical model proposed is found sustainable to predict a quantitative extent of tumor selective accumulation of medicinal nanoparticles in human retinoblastoma cells as long as these NPs are amphiphilic and membranotropic agents with a marked mode for permeability into the target cell.

The C60-fullerene based Me2+ — carrying-n-releasing members of PMC16 family fit the above specified requirements. So our RB/RT proliferation "rate gap" focused computational technique might make a difference in optimization of the preclinical research program for these and related pharmacophores.

References

1. Buchachenko AL. Magneto-Biology and Medicine. New York: Nova Biomedical Publ., 2015.

2. Bukhvostov AA, Dvornikov AS, Ermakov KV, Kurapov PB, Kuznetsov DA. Retinoblastoma: magnetic isotope effects might make a differece in the current anti-cancer research strategy. Acta Medica (Prague). 2017; 60 (2): 93-6.

3. Bukhvostov AA, Dvornikov AS, Ermakov KV, Kuznetsov DA. Retinoblastoma case: shall we get a paramagnetic trend in chemotherapy? Arch Cancer Res. 2017; 5 (4): 158-62.

4. Bukhvostov AA, Pavlov KA, Ermakov KV Sidoruk KN, Rybakova IV, Kuznetsov DA, Roumiantsev SA. An atypical ß-like DNA Polymerase of retinoblastoma cells as a target for spin-selective inhibitory cytostatics. J Fund Med Biol (Russian). 2018; 7 (2): 50-3.

5. Orlova MA, Osipova EY, Roumiantsev SA. Effect of 67Zn-nanoparticles on leukemic cells and normal lymphocytes. Brit J Med Med Res. 2012; 2 (1): 21-30.

6. Orlov AP, Orlova MA, Trofimova TP, Kalmykov SN, Kuznetsov DA. The role of zinc and its compounds in leukemia. J Biol Inorg Chem. 2018; 23 (3): 347-62.

7. Siccardi M, Owen A. Towards a computational prediction of nanoparticle pharmakokinetics and distribution. J In Silico & In Vitro Pharmacol. 2016; 2 (1): 8-11.

8. Augsburger JJ, Chow CML, Dyer V, Roussel MF. Translating science into survival. In: Cajjar A, Pappo A, editors. St. Jude's Children's Research Hospital Report. St. Judes CRH Publ.: Memphis, TN — Cincinnati, OH. 2016; 6-55.

9. Zhang J, Benavente CA, McEvoy J, Flores-Otero J, Ding L, Chen X, Ulyanov A. A novel retinoblastoma therapy from genomic and epigenetic analyses. Nature. 2012; 481 (7381): 329-34.

10. Bozic I, Nowak MA. Resisting resistance. Ann Rev Cancer Biol. 2017; 1 (1): 203-21.

11. Bassukas ID. Comparative Gompertzian analysis of alterations of tumor growth patterns. Cancer Res. 1994; 54 (16): 4385-92.

12. Byrne H, Prezidosi L. Modelling solid tumor growth using the theory of mixtures. Math Med Biol. 2003; 20 (4): 341-66.

13. Komarova NL. Mathematical modelling of tumorigenesis: mission possible. Curr Opinion Oncol. 2005; 17 (1): 39-43.

14. Trapp S, Horobin RW. A predictive model for the selective accumulation of chemicals in tumor cells. Eur Biophys J. 2005; 34 (7): 959-66.

15. Udvardi L, Lakatos J, Loewenhaupt RK. Dividing Cell. In Vitro-In Silico Models. Szeged, Budapest: Alba Regia, 2017.

16. Lehman RJ, Waugh TS, Rattenau KR, Bielka H. An expanding compartment mode to implement a guest probe diffusion input adopted by the comprehensive Gunault muxtures theory. In: Sieliwanowicz B, Martell SJ, edsitors.Combinational Dynamics in Systems Theory. Sydney—Melbourne—Perth: Adler & Adler Publ., 2018; 116-37.

17. Lamprecht A, Pellecker J. Anti-cancer nano-size agents: targeting paths and pharmacokinetics. In: Lamprecht A, editor. Nanotherapeutics. Drug Delivery Concepts in Nanoscience. NY—London—Singapore: Pan Stanford Publ., 2018; 92-101.

18. Altinok A, Gonze D, Levi F, Goldbeter A. An Automaton model for the cell cycle. Interface Focus. 2011; 1 (1): 36-47.

19. Delbreaux J, Pitot CA. Mathematical model in a new drug preclinical trial. Predictive power and limitations. Leuven— Ghent—Antwerp: Leuven University Press, 2018.

Литература