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 ²⁵Mg²⁺, ⁴³Ca²⁺, ⁶⁰Co²⁺, ⁶⁷Zn²⁺ 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 (t – τ)[1 – D(t)] ,

where x(t) is population density at time t, λ is the cell birth rate, τ 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 (Tm) 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_m/T_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 τ > δ > τ/2:

T = τω/|τ – (δ + ω)| .

Here, δ is the duration of the period in which the drug effective, and ω is the period in which the drug dosage is below efficiency. The point, τ = δ + ω, 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 phenomenon) and a routine pharmacokinetic key points prediction.

Inter-specimen covariant models

A. Pharmacokinetic model

form. 1

B. Non-vectoral covariative model

CL_j = [θ₃ • OCC1 + θ₄(WT – 75)] • exp (η_CLj)

V_cj = [θ₁ – (GFR – 80) • θ₂] • exp (η_Vj)

k_12j = (θ₅) • exp (ηk_12j)

C. Population dynamics model

form. 2

D. Population parameters

θ₁ = 19,5

θ₂ = 0,198

Gompertzian model

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 τ_φ have been discretized in several elementary age intervals α ∈ {1,…, N_ϕ } where N_ϕ is an integer such as τϕ = dt • N_ϕ. Here dt is the time step of the cell cycle model. The cell density n_α,ϕ at age α in phase ϕ is governed by:

form. 3

In this equation, ϕ ∈ {G₁, S, G₂ M, G₀, Apoptosis and α ∈ {1,…,Nϕ}.

P_α,ϕ 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 G₁ phase.

Furthermore, noting that Σ_α,ϕ n_α,ϕ is constant, so we can sum to obtain an expression for the pressure field of the form:

form. 4

The computer program starts from an initial distribution of neighboring RB and RT cells in each state {α, ϕ}. 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:

form. 5

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 (G₁, S, G₂ and M), this might be re-formalized as:

form. 6

along with the initial conditions (ni = 0 )_1≤i≤I .

Cell death rates in phases are noted d_ϕ and transition rates between phases, assumed to be time-periodic. K_i→i+1. Phase i (1 ≤ i ≤ I) may be one of the classical four G₁, S, G₂ and M, but also an aggregated phase such as S-G₂, or even a single proliferating phase G₁-S-G₂-M, or, on the contrary, a subdivision inside a phase, e.g., pre- or post-restriction point in G₁; the equation describes the evolution of the densities n_i (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 ʄ for cytotoxicity:

form. 7

whereas λ, ν, ɛ_D, α, B_max, H, φ_B, ϒ_B, D₅₀ are positive constants, identified on tumor growth curves or from literature data [8, 15, 18], or else estimated.

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

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

form. 8

where

form. 9

and λ, μ, ɛ_C, α, β, ϒ, Z_eq, F, φ_A, ϒ_A, C₅₀ 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) (Ф representing a “tap on-tap off” function), and health tissue homeostasis by a linear system showing a stable focus at Z_eq, A_eq = β^-1 ( ϒ – αZ_eq ), perturbed by the drug cytotoxicity function which comes to strengthen the natural self-regulation coefficient α.

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

form. 10

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:

form. 11

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 is ʄ, rate of escape at mitosis towards the siding phase Q, ʄ 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 Me²⁺- 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).

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:

A_C = K_d [tgα_RB / (tgα_RB – tgα_RT) ] ,

where K_d 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 C₆₀-fullerene based Me²⁺ — 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.