1 Chlamydia Trachomatis

1.1 Development Cycle

Chlamydia trachomatis is a major cause of genital and eye infections in humans [2, 7]. There were an estimated 131 million new cases of C. trachomatis genital infection in 2012, making it the most common cause of sexually transmitted bacterial infection [14]. This pathogen also causes trachoma, which is the world’s leading cause of preventable blindness. Twenty one million people are estimated to have active trachoma [17] of which 1.9 million individuals are blind or have severe visual impairment [22]. Trachoma has been targeted for elimination by the World Health Organization.

The C. trachomatis replicates by means of an unusual 2-day developmental cycle that takes place within a human, or other eukaryotic, host cell (Fig. 1) [1, 11, 13]. Unlike most bacteria, Chlamydia exists in two distinct morphologic forms, each with specialized function. The infectious form, called an elementary body (EB), does not proliferate but infects a host cell by entering it from extracellular space and immediately changing into a reticulate body (RB), the replicating form of the bacterium to be enclosed by a membrane-bound compartment known as the chlamydial inclusion. Within the inclusion, the RB unit divides repeatedly by binary fission to expand the RB population. After an initial period of no conversion (known as a conversion holiday), an individual RB may convert into an EB. This conversion event occurs asynchronously, so that some RBs are converted into EBs, while others continue to divide. At later times of the life cycle, an inclusion would contain a mixture of RBs and EBs of different vintages. This mixture of chlamydial forms is released into the extracellular space when the host cell lyses [1, 10]. The RB-to-EB conversion is a critical step in the developmental cycle since only EBs can survive outside the host cell to infect new host cells to further spread the infection. The cartoon of Figure 6.2 in [20] summarizes this two-form characterization of the C. trachomatis life cycle.

Until recently, these developmental events of intracellular Chlamydia infection had not been quantified. Conventional electron microscopy of a Chlamydia-infected cell analyzes a (planar) section through the inclusion, and thus reveals the relative proportions, but not accurate numbers, for each developmental form [3, 9, 16]. Using a novel three-dimensional electron microscopy approach known as serial block-face scanning electron microscopy (SBEM) [5, 6, 12], the labs of Tan and Sütterlin recently performed a comprehensive analysis of chlamydial inclusions. The analysis quantified the changing number of RBs and EBs in the inclusions over the course of the developmental cycle [11]. Their study also provided extensive quantitative data on RB replication and RB-to-EB conversion because it was able to identify and quantify RBs undergoing replication or conversion. Other features observed include the delayed and asynchronous nature of RB-to-EB conversion, and a 6-fold decrease in average RB size as the RB population expanded through repeated divisions by binary fission. Three-dimensional SBEM reconstructions of representative Chlamydia inclusions of infected host cell can be found in Figure 3 of [11].

Taking into account the correlation between small RB size and the onset of RB-to-EB conversion, Lee et al. proposed a size control mechanism in which RBs reduce in size through repeated division, and cannot convert into an EB until they are below a size threshold. Adhering to such a size control conversion mechanism, growth curves of Chlamydia populations can be replicated by a probabilistic model utilizing Gamma distributions for each transition time between (division and conversion) stages [11]. Subsequently, a different theoretical investigation of the same life cycle was also initiated to seek a bio-theoretical basis of the observed life cycle by modeling RB proliferation as a simple birth (proliferation) and death (conversion) process. With conversion rate as the control, it has been shown in [8, 21] that the bacteria’s mean growth dynamics and delayed RB-to-EB conversion are consequences of the bacterium’s strategy to maximize the production of infectious progeny and, thus, the spread of the intracellular infection to new host cells.

1.2 Maximum Spread of Infection and Selective Pressure

With a high capacity for the RB-to-EB conversion (relative to the RB proliferation rate), allowing more RBs to divide at an early stage would generate more RBs for later conversion to more infectious EBs to be discharged when the host cell lyses. With its initial conversion holiday, the Chlamydia development cycle would increase the terminal EB population for a larger spread of the bacterial infection. As the expected bacterial response to the pressure of species survival in the biological world of competitive exclusion, maximizing the terminal EB population should involve some form of delayed conversion.

To assess the validity of this posit, two-form deterministic models for the evolution of mean RB and EB populations (deduced from the birth and death process model in [8]) is formulated and analyzed in [8, 21] for a typical chlamydial inclusion. With the (per unit RB) conversion rate u(t) as the control, these models seek the optimal u(t) that maximizes the EB population at the terminal time T when the host cell lyses. For the biologically relevant conversion capacity range, the optimal maximizing strategy does in fact start with a conversion holiday. The small size threshold for conversion eligibility proposed in [11], not stipulated in the models of [8, 21] is taken simply as an instrument for implementing this initial period of no conversion and plays no role in arriving at the optimal conversion strategy.

In reality, the lysis time of an infected cell varies. Data collected in [11] and elsewhere also show variability in other developmental features. For this reason, a probabilistic model with an uncertain lysis time has also been investigated (with the error in [20] corrected in [8]), again with the size threshold eligibility for conversion playing no role in a similar bang-bang control strategy for maximizing the terminal EB population.

1.3 RB Growth after the Conversion Holiday

To minimize the complexity of analysis, optimal control models were first formulated for two (lumped) forms (RB and EB) of the Chlamydia population and hence without features such as repeated RB divisions or a small size threshold conversion eligibility. It is therefore not surprising that the prediction from such models may not reflect all aspects of the life cycle shown in available empirical data. That they all predict a conversion holiday at the start of the life cycle (so that the observed conversion can be interpreted as the bacterium’s optimal strategy for species’ competitive survival) is already a pleasant surprise.

Still, it has been noted in [8] a qualitative difference between predictions by our optimal control models and the available data after the initiation of the RB-to-EB conversion. More specifically, data collected in [11] show a continual growth of the RB population after the initiation of the RB-to-EB conversion at least for a period. In contrast, the optimal control models of [8, 21] all show an immediate and continual decline of the (lumped) RB population after the onset of the RB-to-EB conversion. It is the goal of the present effort to show how this deficiency can be remedied by a multi-division model with the RB-to-EB conversion possible only when threshold size eligibility is met. To set the present approach in proper context, the basic two-form optimal control model for the time evolution of the two (lumped) Chlamydia populations RB and EB and the relevant results deduced from the model will be summarized in the next section.

2 The Basic Optimal Control Model

2.1 The Optimal Control Problem

For this simplest two-form model of [21], the growth of the (mean) RB population R(t) and the EB population E(t) in physical time t is determined by the two IVP

$$\begin{aligned} R^{\prime } & = (\alpha -u)R,\ \ \ \ R(0)=N, \end{aligned}$$
(1)
$$\begin{aligned} E^{\prime } & = uR,\ \ \ \ E(0)=0, \end{aligned}$$
(2)

respectively, where \(\alpha\) is the natural RB growth rate, u(t) the RB-to-EB conversion rate, and \((\ \ )^{\prime }=\frac{\mathrm{{d}}(\ \ )}{\mathrm{{d}}t}\). For R(t) and E(t) to be continuous and piecewise differentiable, we limit admissible u(t) to the set \(\varOmega\) of piecewise smooth (PWS) functions. The bacterium’s strategy for competitive survival is to maximize the EB population at the terminal time T when the host cell lyses. With the conversion rate u(t) as the control at its disposal, the optimization problem is to choose u(t) to maximize \(E_{T}=E(T)\):

$$\begin{aligned} \max _{\,u \,\in \varOmega }\ \bigg [\,E_{T}=\int _{0}^{T}u(t)R(t)\,\mathrm{{d}}t \bigg ], \end{aligned}$$
(3)

subject to the growth dynamics and the initial condition (1) for R(t). The terminal time T is taken to be a known cell lysis age, though it may also be determined alternatively by the host cell’s capacity constraint

$$\begin{aligned} E(T)+m_{\text{c}}\,R(T)=P_{\text{c}} \end{aligned}$$
(4)

for some prescribed capacity \(P_{\text{c}}\) with \(R_{T}=R(T)\). Both cases have been previously analyzed in [8, 21]. For the intended purpose of the present investigation, it suffices to consider the case of a known cell lysis time so that the terminal time T is a prescribed quantity.

Given the reality of a limit to the capacity to convert, there is an upper bound on the one-way RB-to-EB conversion. The optimal choice \(u_\mathrm{{op}}(t)\) among all u(t) in the admissible set \(\varOmega\) is required to satisfy the inequality constraints

$$\begin{aligned} 0\leqslant u(t)\leqslant u_{\max }. \end{aligned}$$
(5)

We will focus on the biologically more interesting case

$$\begin{aligned} u_{\alpha }\equiv u_{\max }-\alpha >0 \end{aligned}$$
(6)

in the subsequent development. Otherwise, the optimal control would necessarily be \(u_{\max }\) from the start since the RB population would continue to grow even with the maximum RB-to-EB conversion.

2.2 The Optimal Conversion Strategy

The appropriate method for finding the optimal conversion rate \(u_\mathrm{{op}}(t)\) that maximizes the (mean) terminal EB population \(E_{T}\) is the maximum principle [4, 15, 18]. By this method, it is shown in [8, 21] that for the biologically realistic case of \(u_{\alpha }>0\) the optimal conversion strategy is

Proposition 1

$$\begin{aligned} u_\mathrm{{op}}(t)=\left\{ \begin{array}{c} 0 , \ \ \ \ \ \ \ \ \ \ \ 0<t<t_{\text s}, \\ u_{\max , }\ \ \ \ \ \ \ t_{\text s}<t\leqslant T , \end{array} \right. \end{aligned}$$
(7)

where the switch point \( t_{\text s}\) is given by

$$\begin{aligned} {t_{\text s}} =T-\frac{1}{u_{\alpha }}\ln \left( \frac{u_{\max }}{\alpha }\right) <T \end{aligned}$$
(8)

given \(0<\alpha <u_{\max }\) with \( t_{\text s}\downarrow 0\) as \(u_{\max }\downarrow \alpha\).

Remark 1

Note that the switch point occurs prior to the terminal time.

A conversion strategy of the form (7) involving a switch from one extreme conversion rate to the other extreme rate with a finite jump discontinuity at the instant \( t_{\text s}\) is known as a bang-bang control in the control literature [4, 18].

2.3 Terminal Chlamydia Populations

With the optimal bang-bang control (7), it is straightforward to use the growth dynamics (1) and (2) to determine the corresponding R(t) and E(t) to be

$$\begin{aligned} R_\mathrm{{op}}(t)=\left\{ \begin{array}{c} N \text{e}^{\alpha t},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0\leqslant t\leqslant t_{\text s}, \\ N\text{e}^{\alpha t_{\text s}}\text{e}^{-u_{\alpha }(t-t_{\text s})},\ \ \ \ \ t_{\text s}\leqslant t\leqslant T, \end{array} \right. \end{aligned}$$

and

$$ \begin{aligned} E_\mathrm{{op}}(t)=\left\{ \begin{array}{c} 0 ,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0\leqslant t\leqslant t_{\text s,} \\ \frac{u_{\max }}{u_{\alpha }}N {\text {e}}^{\alpha t_{\text s}}\left[ 1- \text e^{-u_{\alpha }(t-t_{\text s})}\right] , \ \ \ \ \ t_{\text s}\leqslant t\leqslant T \end{array} \right. \end{aligned} $$
(9)

with the optimal time \( t_{\text s}\) for the onset of the RB-to-EB conversion given by (18) in terms of the known cell lysis time T.

The optimality of (9) has been discussed in [8, 21]. For instance, it is shown there that \(E_\mathrm{{op}}(T)\) is generally large compared to \(\left[ E_{T}\right] _{u(t)=\alpha }\).

2.4 The RB Populations After Conversion Holiday

That the optimal conversion strategy from our optimal control model for maximizing the terminal EB population turns out to be a bang-bang control (7) for \(\alpha <u_{\max }\) is gratifying. Without any presumption of an initial conversion holiday or any enabling mechanism (such as a size-threshold for conversion) to promote its occurrence, the optimal conversion strategy deduced unequivocally anticipates the experimental findings in [11]. If nothing else, they help to buttress our expectation that species survival is at work to result in the observed developmental cycle.

The relevant data on the two (lumped) populations (averaged over inclusions) summarized in Table 1 have been excerpted from the pie chart of Figure 1 in [11]. The data support qualitatively the theoretical optimal conversion strategy (7). The pie-charts in Figure 1 of [11] show only RB in the infected cells 20 hours post infection (hpi). With data collected every 4 hpi, the conversion to EB begins sometimes around 24 hpi [5]. This is well after the infecting EB entering the host cell, over a time period nearly half way through the development cycle.

Table 1 Mean count of RB and EB (averaged over inclusions)
Full size table

The total EB population rises quickly to 192 units at 28 hpi indicating a maximum possible conversion rate that is much greater than the natural RB proliferation rate. They also show a sharp rise of EBs suggesting a nearly discontinuous conversion rate around \( t_{\text s}\lesssim 24\) hpi. Furthermore, with the total RB population declining sharply after reaching a peak of 507 units at 32 hpi, the corresponding gain in EB units becomes slower in the data collected (up to 40 hpi). These observations provide further support for conversion being proportional to the RB population size as modeled by (2) and with a conversion rate upper bound \(u_{\max }R\) considerably higher than the natural growth rate \(\alpha R\) for the RB form.

It would be natural at this point to make use of parameter values (for \(\alpha\), \(u_{\max }\), etc.) estimated in [11] to calculate relevant quantities such as the switch point \(t_{\text s}\), terminal time T, and other information pertaining to the mean populations of interest. However, the available data from Table 1 excerpted from [11] suggest that the basic model of this section (or any other two-form models of [8, 21]) still needs further refinement before embarking on this task. From the growth dynamics of the RB population modeled by (1), we see for the relevant range \(\alpha <u_{\max }\) of interest here that the RB population grows exponentially prior to the onset of conversion (since \(u_\mathrm{{op}}(t)=0\) during the conversion holiday at the start) but begins to decline immediately after the switch point \( t_{\text s}\) (\(\simeq 24\) hpi) since \(R^{\prime }=\alpha -u_{\max }<0\) (for \(t> t_{\text s}\)). On the other hand, the total RB population reported in Table 1 continues to increase for a period after the appearance of EB particles around 24 hpi. That total only starts to decline after reaching a maximum of 507 units at 32 psi. This important qualitative difference is likely due to the lumped form simplification in the basic model (and other two-form models). Without multiple divisions to reach the conversion eligibility threshold size, two-form models are not sufficiently fine-grained for matching the reported data beyond the conversion holiday phase of the bacterial growth. In this article, we investigate a multi-stage proliferation model with the size eligibility for conversion for the possibility of rectifying the discrepancy in the RB growth during the conversion phase of the life cycle.

3 Multi-step Growth and Single-Step Conversion

3.1 A Multi-stage Proliferation Model

Though artificially lumped into a single RB population group (for further proliferation or conversion), intermediate RB divisions have been identified and reported separately in [11] for their distinctive appearance on electron micrographs with RB-to-EB conversion possible only after n (seemingly less than 6) successive divisions. A more realistic model for C. trachomatis proliferation and differentiation would work with an n-form RB model for \(\left\{ R_{k}(t)\right\}\), \(k=1,2,\cdots ,n\) (as in the schematic diagram of Fig. 5a of [11]) with \(R_{j}(t)\) being the population of the RB having had (\(j-1\))th divisions and known as the vintage j RB subpopulation. At any instant of time in this model, some RB units of the vintage k may divide by binary fission to results in two smaller size RB of the \(R_{k+1}\) vintage. At stage n, the RB of the type \(R_{n}\) would be of sufficiently small size to be eligible for conversion into EB units. The growth rates of the n groups of RB populations are summarized mathematically by the following differential equations:

$$\left\{ \begin{aligned} R_{1}^{\prime } & = -\,\alpha _{1}R_{1}\equiv f_{1}(\mathbf {R,\alpha }),\ \ \ \ \ \nonumber \\ R_{2}^{\prime } & = 2\alpha _{1}R_{1}-\alpha _{2}R_{2}\equiv f_{2}(\mathbf {R},{\varvec{{\alpha }}}),\ \ \nonumber \\& \cdots\ \nonumber \\ \ R_{n}^{\prime } & = 2\alpha _{n-1}R_{n-1}-uR_{n}+\alpha _{n}R_{n}\equiv f_{n}(\mathbf {R,\alpha }),\ \ \ \ \end{aligned}\right.$$
(10)

where \({\mathbf {R}}=\left( R_{1},\cdots,R_{n}\right)\) and \(\ \ \ \mathbf {\alpha }=\left( \alpha _{1},\cdots,\alpha _{n}\right)\), starting with some initial populations

$$\begin{aligned} R_{1}(0)=N,\ \ \ \ \ R_{k}(0)=0 , \ \ \ \ \ 1<k\leqslant n. \end{aligned}$$
(11)

When repeated binary fissions reduce the RB form to a conversion eligible size \(R_{n}\), our n-form model characterizes its further development by a growth rate \(\alpha _{n}R_{n}\) for brevity. We also omitted for a proof of concept investigation the possibility of any “vintage k” RB units growing in size to become a member of vintage \(\ k-1\) RB subpopulation (allowed in Fig. 5a of [11]).

3.2 Maximization of Terminal EB Population

The accumulation of EB from conversion at stage n is characterized by

$$\begin{aligned} E^{\prime }=uR_{n},\ \ \ \ \ \ \ E(0)=0 \end{aligned}$$
(12)

with the finite conversion capacity characterized by the bounds in (5) on the per unit RB conversion rate u(t). The total accumulation of EB units at the prescribed terminal time T is then given by

$$\begin{aligned} E_{T}\equiv E(T)=\int _{0}^{T}u(t)R_{n}(t)\,\mathrm{{d}}t \end{aligned} . $$
(13)

Subject to growth dynamics (10)–(12) and the conversion rate constraints (5), the optimal conversion strategy is to choose \(u(t)=u_\mathrm{{op}}(t)\) to maximize (13)

$$\begin{aligned} \max _{0\,\leqslant \,u\,\leqslant \,u_{\max }}[E_{T}]=\left[ E_{T}\right] _{u=u_\mathrm{{op}}(t)}. \end{aligned}$$
(14)

3.3 The Maximum Principle

A method of solution for the optimal control problem described in the previous subsection is by the maximum principle [4, 15, 18, 19]. For this approach, the relevant Hamiltonian function is

$$\begin{aligned} H[u] & = uR_{n}+\sum _{k=1}^{n}\lambda _{k}f_{k}(\mathbf {R,\alpha }) \nonumber \\ & = uR_{n}+\lambda _{n}\left ( 2\alpha _{n-1}R_{n-1}-\left( u-\alpha _{n}\right) R_{n}\right) +\cdots +\,\lambda _{1}\left( -\alpha _{1}R_{1}\right) . \end{aligned}$$
(15)

From this we obtain the adjoint ODEs and terminal conditions for the adjoint variables \(\left\{ \lambda _{k}(t)\right\}\):

$$\begin{aligned} \lambda _{k}^{\prime }=-\frac{\partial H}{\partial R_{k}},\ \ \ \ \ \ \ \lambda _{k}(T)=0, \ \ \ \ \ \ k=1, \cdots,n. \end{aligned}$$
(16)

The optimal conversion rate \(u_\mathrm{{op}}(t)\) is then chosen by the condition of optimality

$$\begin{aligned} \max _{0\leqslant u\leqslant u_{\max }}\left[ H[u]\right] =H[u_\mathrm{{op}}]. \end{aligned}$$
(17)

3.4 The Adjoint Problem

Since H[u] is linear in \(R_{k}\), the adjoint ODEs do not depend on the state variables with

$$\begin{aligned} \lambda _{1}^{\prime } & = \alpha _{1}\left( \lambda _{1}-2\lambda _{2}\right) ,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \lambda _{1}(T)=0, \\& \cdots \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ \lambda _{n-1}^{\prime } & = \alpha _{n-2}\left( \lambda _{n-2}-2\lambda _{n-1}\right) ,\ \ \ \ \ \ \ \lambda _{n-1}(T)=0 , \end{aligned}$$

and

$$\begin{aligned} \lambda _{n}^{\prime }=u_{\alpha }\lambda _{n}-u,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \lambda _{n}(T)=0 \end{aligned}$$

with \(u_{\alpha }=u-\alpha _{n}\). Once u(t) is prescribed, the system can be solved recursively backward from \(k=n\), \(k=n-1,\ \cdots \) to \(k=1\).

With \(\lambda _{n}(T)=0\) and \(u_{\max }>\alpha _{n}\) (for \(u_\mathrm{{op}}(t)=u_{\max }\) otherwise), we have \(\lambda _{n}^{\prime }(T)<0\) for any \(u(T)>0\). And, by the continuity of \(\lambda _{n}(t)\), \(\lambda _{n}(t)<1\) for \(t> t_{ \text s}\) for some \( t_{ \text s}\) \(<T\). The condition of the optimality then requires \(u_\mathrm{{op}}(t)=u_{\max }\) for t in \(( t_{ \text s},T]\):

$$\begin{aligned} \max _{0\,\leqslant \,u\leqslant \,u_{\max }}\left[ uR_{n}(1-\lambda _{n})+\left\{ \cdot \cdot \cdot \right\} \right] =H\left[ u_{\max }\right] , \ \ \ \ t_{ \text s} <t\leqslant T \text {.} \end{aligned}$$

In that case, we have

$$\begin{aligned} \lambda _{n}(t)=\frac{u_{\max }}{u_{\alpha }}\left (1-\text e^{-u_{\alpha }(T-t)}\right) , \ \ \ \ \ t_{\text s}<t\leqslant T \end{aligned} $$

with the remaining adjoint functions determined by solving the relevant IVP for \(\left\{ \lambda _{k}(t)\right\}\) recursively, backward starting from \(k=n-1\). These solutions will not be listed here as they will not be needed in the subsequent development.

3.5 The Optimal Conversion Strategy

The choice of \(u_\mathrm{{op}}=u_{\max }\) ceases to be appropriate when \(\lambda _{n}(t)>1\) with the condition \(\lambda _{n}(t_{\text s})=1\) delimiting its range of applicability to \( t_{\text s}<t\leqslant T\). This gives for \(u_{\alpha }=u_{\max }-\alpha _{n}>0\),

$$\begin{aligned} t_{\text s}=T-\frac{1}{u_{\alpha }}\ln \left( \frac{u_{\max }}{\alpha _{n}}\right) <T\text {.} \end{aligned}$$
(18)

With \(u_{\max }\) no longer optimal for \(t< t_{\text s}\), we should consider \(u(t)<u_{\max }\) for \(t< t_{\text s}\). It is not difficult to see that for any constant \({\bar{u}}>0\) (and \(\alpha _{n}>0\)), the choice \(u_\mathrm{{op}}(t)\) \(={\bar{u}}\) would result in \(\lambda _{n}({\bar{t}})>1\) for some time \({\bar{t}}_{s}>0\). This suggests that we consider \(u_\mathrm{{op}}(t)\) \(=0\), the lower bound for u(t) specified in (5), for \(t< t_{\text s}\). For this choice of u(t), we have for the adjoint equation for \(\lambda _{n}\) and the continuity requirement \(\lambda _{n}( t_{\text s})=1\),

$$\begin{aligned} \lambda _{n}(t)=\text e^{\alpha _{n}( t_{\text s}-t)}, \ \ \ \ \ t\leqslant t_{\text s}. \end{aligned}$$
(19)

With \(\lambda _{n}(t)>1\) for \(t< t_{\text s}\), \(u_\mathrm{{op}}(t)\) \(=0\) continues to be optimal in the range of \(0\leqslant\) \(t< t_{\text s}\) by the optimality condition. This result complements the previous result for the range \( t_{s}<t\leqslant T\) to give the following conclusion:

Proposition 2

For the biologically relevant range of \(u_{\max }>\alpha _{n}\), the optimal RB-to-EB conversion strategy is the bang-bang control

$$\begin{aligned} u_\mathrm{{op}}(t)=\left\{ \begin{array}{c} 0, \ \ \ \ \ \ \ \ \ 0\leqslant t<t_{\text s}, \\ u_{\max },\ \ \ \ \ t_{\text s}<t\leqslant T \end{array} \right. \end{aligned}$$
(20)

with the switch time \( t_{\text s}\) given by (18).

The remaining adjoint equations for \(\left\{ \lambda _{n-1}(t),\cdots ,\lambda _{1}(t)\right\}\) may then be solved for this optimal conversion rate for all t in [\(0, t_{\text s}\)) with \(u_\mathrm{{op}}(t)=0\) and with the continuity requirements for these functions at the switch point \( t_{\text s}\) providing a terminal condition for each of the adjoint equations.

4 Evolution of the RB Populations (\(n=3\))

4.1 The Switch Point

To examine the consequence of the new multi-step model on the RB population, we work out here the \(n=3\) case. With \(u_{\alpha }=u_{\max }-\alpha _{3}>0\), we have from (18)

$$\begin{aligned} t_{\text s}=T-\frac{1}{u_{\alpha }}\ln \left( \frac{u_{\max }}{\alpha _{3}}\right) <T \end{aligned}, $$
(21)

and correspondingly

$$\begin{aligned} \lambda _{3}(t)=\left\{ \begin{array}{c} \text e^{\alpha _{3}(t_{\text s}-t)}, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0\leqslant t\leqslant t_{\text s}, \\ \frac{u_{\max }}{u_{\alpha }}\left ( 1-\text e^{-u_{\alpha }(T-t)}\right ), \ \ \ \ t_{\text s}\leqslant t\leqslant T.\ \end{array} \right. \end{aligned}$$

The remaining two adjoint equations for \(\lambda _{2}(t)\) and \(\lambda _{1}(t)\) may then be solved for this optimal conversion rate for all t in [0, T] with the continuity requirements for these functions at the switch point \(t_{ \text s}\) providing a terminal condition for each of these adjoint equations on the range \(0\leqslant t\leqslant t_{\text s}\).

4.2 The Pre-conversion Phase

Prior to the onset of the RB-to-EB conversion at \(t_{\text s}\), the optimal conversion strategy (7) is not to convert. With \(u_\mathrm{{op}}(t)=0\), the equations of state in (36) simplifies slightly to

$$\left\{ \begin{aligned} R_{1}^{\prime } & = -\alpha _{1}R_{1},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ R_{1}(0)=N, \nonumber \\ R_{2}^{\prime } & = 2\alpha _{1}R_{1}-\alpha _{2}R_{2},\ \ \ \ \ \ R_{2}(0)=0,\ \nonumber \\ R_{3}^{\prime } & = 2\alpha _{2}R_{2}+\alpha _{3}R_{3},\ \ \ \ \ \ R_{3}(0)=0. \end{aligned}\right. $$
(22)

These linear ODEs with constant coefficients can be solved exactly to give

$$\left\{ \begin{aligned} R_{1}(t) & = N\text e^{-\alpha _{1}t},\ \ \ \ \ R_{2}(t)=\frac{2\alpha _{1}N}{ \alpha _{2}-\alpha _{1}}\left ( \text e^{-\alpha _{1}t}-\text e^{-\alpha _{2}t}\right ) , \nonumber \\ R_{3}(t) & = \frac{4\alpha _{1}\alpha _{2}N}{\left( \alpha _{1}+\alpha _{3}\right) \left( \alpha _{2}+\alpha _{3}\right) }\left ( \text e^{\alpha _{3}t}+ \frac{\alpha _{1}+\alpha _{3}}{\alpha _{2}-\alpha _{1}}\text e^{-\alpha _{2}t}- \frac{\alpha _{2}+\alpha _{3}}{\alpha _{2}-\alpha _{1}}\text e^{-\alpha _{1}t}\right ) \end{aligned} \right.$$
(23)

for \(0\leqslant t<t_{\text s}\). The corresponding solution for the special case \(\alpha _{2}=\alpha _{1}=\alpha\) of interest is

$$\left\{ \begin{aligned} R_{1}(t) & = N\text e^{-\alpha t},\ \ \ \ \ R_{2}(t)=2\alpha Nt \text e^{-\alpha t},\ \ \ \ \nonumber \\ \ R_{3}(t) & = \frac{4\alpha ^{2}N}{\left( \alpha +\alpha _{3}\right) ^{2}} \left\{ \text e^{\alpha _{3}t}-\left[ 1+(\alpha +\alpha _{3})t\right] \text e^{-\alpha t}\right\} \end{aligned} \right.$$
(24)

for \(0\leqslant t<t_{\text s}\).

While the exact solutions above are useful for characterizing the evolution of the different RB subpopulations, we note parenthetically the following simple result for the special case of \(\alpha _{1}=\alpha _{2}=\alpha _{3}\equiv \alpha\):

Proposition 3

For \(\alpha _{1}=\alpha _{2}=\alpha _{3}\equiv \alpha ,\) the total of all RBs are increasing with time for \(t<t_{\text s}\) with

$$\begin{aligned} R_{1}+R_{2}+R_{3}=R(t)=N\text e^{\alpha t}. \end{aligned}$$
(25)

Proof

Form \(R^{\prime }\equiv (R_{1}+R_{2}+R_{3})^{\prime }\) to get that for \(\alpha _{1}=\alpha _{2}=\alpha _{3}\equiv \alpha\), \(u_\mathrm{{op}}(t)=0\), and \(0\leqslant t<t_{\text s}\),

$$\begin{aligned} R^{\prime }=\alpha R,\ \ \ \ \ R(0)=N , \end{aligned}$$

and

$$\begin{aligned} R_{1}+R_{2}+R_{3}=R(t)=N\text e^{\alpha t}. \end{aligned}$$

For the more general case with \(\alpha _{1}=\alpha _{2}\equiv \alpha\) and \(\alpha _{3}<\alpha\), we have for \(t<t_{\text s}\)

$$\begin{aligned} R^{\prime }=\alpha R-(\alpha -\alpha _{3})R_{3}, \end{aligned}$$

and it shows an exponential growth in the total RB population R(t) which is somewhat reduced by the diminutive contribution from \(R_{3}\). The actual exact total RB population is

$$\begin{aligned} \frac{R(t)}{N}=\frac{4\alpha ^{2}}{\left( \alpha +\alpha _{3}\right) ^{2}} \left\{ \text e^{\alpha _{3}t}-\text e^{-\alpha t}\left[ 1+\left( \alpha +\alpha _{3}\right) t\right] \right\} +\text e^{-\alpha t}\left( 1+2\alpha t\right) \end{aligned}$$

which is reduced to \(R(t)=N \text e^{\alpha t}\) for \(\alpha _{3}=\alpha\). A typical set of \(\left\{ R_{k}(t)\right\}\) and their sum R(t) are shown in Fig. 1 for \(\alpha _{1}=\alpha _{2}=1\) and \(\alpha _{3}=0.25\) (corresponding to t measured in units of the common growth rate constant \(\alpha\)).

It would seem therefore that the collection of the different RB populations is likely to be increasing immediately after the onset of the RB-to-EB conversion unless the system has a sufficiently high conversion rate capacity relative to the growth rate \(\alpha _{3}\). If \(u_{\alpha }\) should be so large that \(u_{\alpha }R_{3}\) becomes large compared to \(2\alpha _{2}R_{2}\) immediately after the switch time \(t_{\text s}\), the dominant component \(R_{3}(t)\) as well as R(t) may decrease with time. How this qualitative supposition should play out quantitatively will be seen in the next subsection.

Fig. 1
figure 1

Graphs of \(R_{k}(t)\) and R(t) for \(0\leqslant t\leqslant t_{\text s}\) with \(\ \alpha _{1}=\alpha _{2}=1\) and \(\alpha _{3}=0.25\)

Full size image

4.3 The Conversion Phase of the Life Cycle

For \(t>t_{\text s}\), Proposition 2 requires \(u_\mathrm{{op}}(t)=u_{\max }\). While the exact solution for the corresponding IVP for the state equations

$$\left\{ \begin{aligned} R_{1}^{\prime } & = -\alpha _{1}R_{1},\ \ \ \ \ \ \ R_{1}(0)=N \nonumber, \\ R_{2}^{\prime } & = 2\alpha _{1}R_{1}-\alpha _{2}R_{2},\ \ \ \ \ \ \ \ \ \ R_{2}(0)=0, \\ R_{3}^{\prime } & = 2\alpha _{2}R_{2}-u_{\alpha }R_{3},\ \ \ \ \ \ \ \ R_{3}(0)=0 \end{aligned} \right. $$
(26)

is straightforward, we record here only the solution for the special case \(\alpha _{1}=\alpha _{2}=\alpha\). For this special case, we have

$$\begin{aligned} R_{1}(t) & = r_{1}\text e^{-\alpha (t-t_{\text s})},\ \ \ \ \ R_{2}(t)=\left( r_{2}+2\alpha r_{1}\left( t-t_{\text s}\right) \right) \text e^{-\alpha \left( t-t_{\text s}\right) }, \\ R_{3}(t) & = r_{3}\text e^{-u_{\alpha }(t-t_{\text s})}-\frac{\left( 2\alpha \right) ^{2}r_{1}}{\left( u_{\alpha }-\alpha \right) ^{2}}\left( \text e^{-\alpha (t-t_{\text s})}-\text e^{-u_{\alpha }(t-t_{\text s})}\right) \\& \, \ \ \ +\frac{2\alpha \left[ r_{2}+2\alpha r_{1}(t-t_{\text s})\right] }{u_{\alpha }-\alpha }\text e^{-\alpha (t-t_{\text s})}-\frac{2\alpha r_{2}}{u_{\alpha }-\alpha } \text e^{-u_{\alpha }(t-t_{\text s})},\end{aligned}$$

where we have applied the continuity requirements \(R_{k}(t_{\text s})=r_{k}=\left[ R_{k}(t_{\text s})\right] _{u=0}\) with

$$\begin{aligned} r_{1} & = N \text e^{-\alpha t_{\text s}},\ \ \ \ r_{2}=2\alpha Nt_{\text s}\text e^{-\alpha t_{\text s}},\ \ \ \ \ \end{aligned}$$
(27)
$$\begin{aligned} r_{3} & = \frac{4\alpha ^{2}N}{\left( \alpha +\alpha _{3}\right) ^{2}}\left\{ \text e^{\alpha _{3}t_{\text s}}-\left[ 1+\left( \alpha +\alpha _{3}\right) t_{\text s}\right] \text e^{-\alpha t_{\text s}}\right\}. \end{aligned}$$
(28)

The expression for \(R_{3}(t)\) is rather complex to readily reveal its evolutionary behavior. On the other hand, the rate of change of the sum R(t) at \(t_{\text s}{\small +}\) is more informative. The expression

$$\begin{aligned} R^{\prime }(t_{\text s}{\small +})=\alpha R(t_{\text s})-\left( u_{\alpha }+\alpha \right) R_{3}(t_{\text s})=\alpha \left( r_{1}+r_{2}\right) -u_{\alpha }r_{3} \end{aligned}$$
(29)

shows that \(R^{\prime }(t_{\text s}{\small +})>0\) if

$$\begin{aligned} u_{\alpha }=u_{\max }-\alpha _{3}<\frac{r_{1}+r_{2}}{r_{3}}\equiv U_{3} \end{aligned}$$
(30)

with

$$\begin{aligned} U_{3}=\left( \frac{\alpha +\alpha _{3}}{2}\right) ^{2}\frac{\text e^{-\left( \alpha +\alpha _{3}\right) t_{\text s}}(1+2\alpha t_{\text s})}{1-\text e^{-\left( \alpha +\alpha _{3}\right) t_{\text s}}\left[ 1+\left( \alpha +\alpha _{3}\right) t_{\text s} \right] }. \end{aligned}$$
(31)

We have then the following important result in meeting the stated goal of this investigation.

Proposition 4

With the switch point \(t_{\text s}\) for the onset of the RB-to-EB conversion determined by (18), the totality of all the RB units for the special case of \(\alpha _{1}=\alpha _{2}=\alpha\) increases with time at least for a period beyond \(t_{\text s}\) if and only if

$$\begin{aligned} u_{\alpha }=u_{\max }-\alpha _{3}<\frac{r_{1}+r_{2}}{r_{3}}\equiv U_{3}. \end{aligned}$$
(32)

With no loss in generality, we may take \(\alpha =1\) (corresponding to time being measured in units of \(\alpha\)), the multiplicative factor \(\text e^{-\left( \alpha +\alpha _{3}\right) t_{\text s}}=\text e^{-\left( 1+\alpha _{3}\right) \,t_{\text s}}\) in the condition (30) limits \(u_{\max }-\alpha _{3}\) to a somewhat narrow range of the relevant rate differential for the continual growth of R(t) beyond the switch point \(t_{\text s}\). That it does occur is seen from the graphs in Fig. 2 for

$$\begin{aligned} \alpha _{1}=\alpha _{2}=1,\ \alpha _{3}=0.25,\ \ u_{\max }=1,\ \ \ T=2 \end{aligned} .$$

An increase in \(u_{\max }\) to 1.85 or higher would violate the condition (32) and lead to an immediate decline in the total RB population after the start of the RB-to-EB conversion as shown in Fig. 3 (for \(u_{\max }=2\)). Nevertheless, Proposition 4 validates theoretically that the size threshold eligibility for the RB-to-EB conversion does offer a mechanism for continual increase in the total RB population beyond the switch point as seen from the data collected in [11] while the optimal conversion strategy for maximizing the spread of the infection continues to be bang-bang as clearly supported by the same data. Furthermore, the graph for R(t) in Fig. 4 also shows that its rise after the onset of conversion is not indefinite. For the particular combination of biological growth and conversion rates, the total RB population peaks and begins to decline well before the terminal time qualitatively similar to the data from [11] shown in Table 1.

Fig. 2
figure 2

Graphs of \(R_{k}(t)\) and R(t) for \(t \geqslant t_{k}\) and n = 3 (with \(\ \alpha _{1}=\alpha _{2}=1\) and \(\alpha _{3}=0.25,\) \(u_{\max }=1\) and \(T=2\)), U3 = 26.587 \(\cdots > u_{\max } - \alpha_{3}\)

Full size image
Fig. 3
figure 3

Graphs of \(R_{k}(t)\) and R(t) for \(t \geqslant t_{k}\) and n = 3 (with \(\ \alpha _{1}=\alpha _{2}=1\) and \(\alpha _{3}=0.25\), \(u_{\max }=2\) and \(T=2\)), U3 = 1.377 78 \(\cdots <u_{\max } - \alpha_{3}\)

Full size image

5 The General Case

The (proof-of-concept) three-stage proliferation model of the previous section confirms that a small size threshold requirement for RB-to-EB conversion could remedy the inconsistency between the simple two-form models of [8, 21] and the empirical data of [11] on the evolution of the total RB population after the onset of conversion. The analysis for the proof-of-concept model can be easily extended to reach a similar conclusion for the general n-stage proliferation model defined by (10)–(12) and (14).

5.1 Prior to the Onset of Conversion

Prior to the onset of the RB-to-EB conversion at \(t_{\text s}\) (as given by (18) with \(u_{\alpha }=u_{\max }-\alpha _{n}\)), the optimal conversion strategy is not to convert. With \(u_\mathrm{{op}}(t)=0\), the equations of state in (36) can be solved exactly. For the case

$$\begin{aligned} \alpha _{1}=\cdots =\alpha _{n-1}=1,\ \ \ \ \ 0<\alpha _{n}<1 \end{aligned}$$

(corresponding to time measured in units of the growth rate \(\alpha _{k}=\alpha ,k<n\)) of interest to us, the exact solutions for individual RB groups are found (by Mathematica or Maple for accuracy) to be

$$\left\{ \begin{aligned} R_{1}(t) & = N \text e^{-t},\ \ \ \ \ R_{2}(t)=\frac{N}{1!}\left( 2t\right) \text e^{-t},\ \ \ \ R_{3}(t)=\frac{N}{2!}\left( 2t\right) ^{2} \text e^{-t}, \nonumber \\ \ R_{4}(t) & = \frac{N}{3!}\left( 2t\right) ^{3} \text e^{-t},\ \ \ \ \ \cdots \ ,\ \ \ \ \ R_{n-1}(t)=\frac{N}{\left( n-2\right) !}\left( 2t\right) ^{n-2} \text e^{-t}, \nonumber \\ R_{n}(t) & = \frac{2^{n-1}N}{\left( 1+\alpha _{n}\right) ^{n-1}}\left ( \text e^{\alpha _{n}t}- \text e^{-t}\sum _{k=1}^{n-1}\frac{\left( 1+\alpha _{n}\right) ^{k-1}}{\left( k-1\right) !}t^{k-1}\right ) . \end{aligned} \right.$$
(33)

Evidently, \(R(t)=\sum _{k=1}^{n}R_{k}\) is dominated by the component \(R_{n}(t)\) and hence grows exponentially during the conversion holiday phase \(0\leqslant t\leqslant t_{\text s}\). In fact, for the special case \(\alpha _{n}=\alpha =1\) so that Proposition 3 holds more generally, we have again \(R(t)=N \text e^{t}\).

5.2 The Conversion Phase of the Life Cycle

For \(t>t_{\text s}\), we have from the bang-bang control (7) that conversion is to be at the maximum allowable rate. The RB subpopulations during the conversion phase are to be determined by the equations of state (36) with \(u_\mathrm{{op}}(t)=u_{\max }\) subject to the continuity requirements

$$\begin{aligned} R_{k}(t_{\text s})=\left[ R_{k}(t_{\text s})\right] _{u=0}=r_{k},\ \ \ \ \ k=1,2,\ \cdots,\ n \end{aligned}$$

with \(\left[ R_{k}(t)\right] _{u=0}\) given in (33).

Evidently, the exact solutions for \(\left\{ R_{k}(t)\right\}\) for the conversion phase of the development are too complex to readily reveal the evolutionary behavior of the total R(t) of all RB subgroups. However, our experience with the \(n=3\) case shows that the rate of change of the sum R(t) at the switch point \(t_{\text s}\) should provide information on how R(t) would evolve shortly after the onset of the RB-to-EB conversion. With

$$\begin{aligned} R^{\prime }(t_{\text s}{\small +})=R(t_{\text s})-u_{\max }R_{n}(t_{\text s})=\sum _{k=1}^{n-1}r_{k}-u_{\alpha }r_{n}, \end{aligned}$$

the expression for \(R^{\prime }(t_{\text s}{\small +})\) shows that the total of all RB units would increase at the start of the RB-to-EB conversion phase if and only if

$$\begin{aligned} u_{\alpha }=u_{\max }-\alpha _{n}<\frac{1}{r_{n}}\sum _{k=1}^{n-1}r_{k}\equiv U_{n} \end{aligned},$$
(34)

where

$$\begin{aligned} U_{n}=\frac{\left[ {\frac{1}{2}} \left( 1+\alpha _{n}\right) \right] ^{n-1} \text e^{-\left( 1+\alpha _{n}\right) t_{\text s}}\sum _{k=1}^{n-1}\left( 2t_{\text s}\right) ^{k-1}/\left( k-1\right) !}{ 1- \text e^{-(1+\alpha _{n})t_{\text s}}\sum _{k=1}^{n-1}\left [ t_{\text s}\left( 1+\alpha _{k}\right) \right ] ^{k-1}/\left( k-1\right) !} \end{aligned} .$$
(35)

We have then the following important result in meeting the state goal of this investigation for the general n-stage model analogous to Proposition 4:

Proposition 5

After the onset of the RB-to-EB conversion at the switch point \(t_{\text s}\), the totality of all RB populations for the special case of \(\left\{ \alpha _{j}=1,j=1, \cdots,n-1,\alpha _{n}<1\right\}\) increases at least for a period beyond \(t_{\text s}\) if and only if (34) is satisfied.

As in the \(n=3\) case, the factor \(\text e^{-\left( 1+\alpha _{n}\right) t_{\text s}}\) in the condition (34) again delimits the range of \(u_{\max }-\alpha _{n}\) for the continual growth of R(t) beyond the switch point \(t_{\text s}\). The important observation here is that

Corollary 1

For a fixed \(t_{\text s}\), the permissible range of the net conversion rate \(u_{\alpha }\) increases with n.

Proof

The corollary follows from the increase of the numerator and the decrease of the denominator as n increases.

For the same \(T=2\) example of Fig. 2, it is seen from Fig. 4 for \(n=4\) that R(t) continues to grow beyond \(t_{\text s}\) (similar to the data from [11] reported Table 1) but now for \(u_{\max }=2.5\) not permissible by the \(n=3\) case. On the other hand, R(t) declines immediately with the onset of conversion for the larger \(u_{\max }=3\) shown in Fig. 5 as it violates the upper bound \(U_{4}\) for the corresponding host cell system. Nevertheless, the result of Corollary 1 further expands the range of \(u_{\max }\) with increasing n for conversion phase RB growth. As such, it further strengthens the adequacy of the n-stage model for rectifying the glaring qualitative disparity of two-form models on the conversion phase growth and declines of the total of all RB subpopulations.

Fig. 4
figure 4

Graphs of \(R_{k}(t)\) and R(t) for \(t \geqslant t_{k}\) and n = 4 (with \(\ \alpha _{1}=\alpha _{2}=1\) and \(\alpha _{3}=0.25\), \(u_{\max }=2.5\) and \(T=2\)), U4 = 2.799 5 \(\cdots >u_{\max } - \alpha_{4}\)

Full size image
Fig. 5
figure 5

Graphs of \(R_{k}(t)\) and R(t) for \(t \geqslant t_{k}\) and n = 4 (with \(\ \alpha _{1}=\alpha _{2}=1\) and \(\alpha _{3}=0.25\), \(u_{\max }=3\) and \(T=2\)), U4 = 2.177 \(\cdots <u_{\max } - \alpha_{4}\)

Full size image

6 A Four-Form Population Model

6.1 Introduction

Since an appropriate n-stage model of the previous section is capable of reproducing the two principal features of the C. trachomatis life cycle, it would be natural to estimate the various growth and conversion rate constants from empirical data to validate the model quantitatively. However, fine-grained information on multiple forms of RB has not been recorded for this purpose. In addition to the initial RB population and final EB population, only RBs in the midst of division (known as dividing RB and denoted by DB) and a new intermediate form IB (called intermediate body) resulting from RB conversion have been recorded separately in [11] and shown in Table 2.

Table 2 Mean counts of Chlamydia subpopulations (averaged over inclusions)
Full size table

A four-form model (of RB, DB, IB, and EB) with a different closure for the two RB subpopulations (RB and DB) has been suggested in [8] for the possible increase of the total RB populations after the onset of the RB-to-EB conversion. An investigation of such a model is to be a separate project given the anticipated tedium (but otherwise routine nature) of the required calculations for the solution of the optimal control problem. We describe briefly in this section the tedious iterative nature of a part of the solution process and the possible transformation of the adjoint problem into one that can be solved by a single pass algorithm.

At any instant of time in this four-form model of [8], some (eligible) RB units of the population R(t) may convert into an Intermediate form I(t) at the rate uR (that would eventually become EB) while a fraction of the remaining R(t) would transition into the DB form at the rate \(\alpha _{1}R(t)\). Each unit of the DB population D(t) is capable of binary division into two new RBs at the rate \(\alpha _{2}D.\) The population I(t) of Intermediate form converts to EB at a rate bI. For the investigation here, it is not necessary to consider the intermediate conversion to IB before arriving at EB. In that case, the growth rates of the C. trachomatis subpopulations are summarized mathematically by the following differential equations:

$$\begin{aligned} R^{\prime } & = -\,(\alpha _{1}+u)R+2\alpha _{2}D,\ \ \ \ D^{\prime }=\alpha _{1}R-\alpha _{2}D, \end{aligned}$$
(36)
$$\begin{aligned} \ \ \ E^{\prime } & = uR \end{aligned}$$
(37)

starting with the initial conditions

$$\begin{aligned} R(0)=N,\ \ \ \ \ D(0)=0,\ \ \ \ \ E(0)=0. \end{aligned}$$
(38)

The problem of determining the PWS optimal RB-to-EB conversion strategy u(t) (to result in the largest possible terminal EB population E(T) at the time T when the host cell lyses) again takes the form

$$\begin{aligned} \max _{0\ \leqslant \ u\ \leqslant \ u_{\max }}\left[ \ E(T)=\int _{0}^{T}uR(t)\mathrm{{d}}t\right] , \end{aligned}$$
(39)

subject to the growth dynamics of the RB and DB of (36) and their associated initial conditions in (11) along with the inequality constraints (5).

6.2 The Optimal Control

6.2.1 Hamiltonian and Optimality

Similar to previous models, the Hamiltonian for the application of the maximum principle is

$$\begin{aligned} H[u] & = uR+\lambda _{1}\left[ -\,(\alpha _{1}+u)R+2\alpha _{2}D\right] +\lambda _{2}\left[ \alpha _{1}R-\alpha _{2}D\right] \\ & = uR\left( 1-\lambda _{1}\right) +\left\{ \cdot \cdot \cdot \right\}, \end{aligned}$$

where terms in \(\left\{ \cdot \cdot \cdot \right\},\) do not involve the control u(t). With the optimality condition (17) and the adjoint boundary condition \(\lambda _{1}(T)=0\), we again arrive at

$$\begin{aligned} u_\mathrm{{op}}(t)=u_{\max }, \ \ \ \ \ t_{\text s}<t\leqslant T \end{aligned}$$

when the adjoint function \(\lambda _{1}(t)\) remains \(<1\). The switch point \(t_{\text s}\) is specified by \(\Lambda _{1}(t_{\text s})\equiv \left[ \lambda _{1}(t_{\text s}) \right] _{u=u_{\max }}=1\). But unlike all previous optimal control models for the life cycle, this condition no longer results in the simple expression (18).

6.2.2 The Adjoint Problem

Adjacent to the terminal time T, the adjoint equations and adjoint boundary conditions are given by

$$\begin{aligned} \lambda _{1}^{\prime } & = \left( \alpha _{1}+u_{\max }\right) \lambda _{1}-\alpha _{1}\lambda _{2}-u_{\max },\ \ \ \ \ \lambda _{1}(T)=0, \end{aligned}$$
(40)
$$\begin{aligned} \lambda _{2}^{\prime } & = \alpha _{2}\left( \lambda _{2}-2\lambda _{1}\right) ,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \lambda _{2}(T)=0 \end{aligned}$$
(41)

for \(t_{\text s}<t\leqslant T\). Unlike the corresponding problems of previous sections, the adjoint function \(\lambda _{1}(t)\) needed for determining \(t_{\text s}\) can no longer be determined separately as in (18). The system (40)–(41) must now be solved simultaneously to obtain

$$\begin{aligned} \lambda _{1}(t)=\Lambda _{1}(t),\ \ \ \ \ \lambda _{2}(t)=\Lambda _{2}(t) \end{aligned}$$

with

$$\begin{aligned} \Lambda _{1}(t) & = \left [ c_{1} \text e^{\mu _{1}(t-T)}+c_{2} \text e^{\mu _{2}(t-T)}\right ] +{\bar{\Lambda }}_{1}, \\ \Lambda _{2}(t) & = \left [ c_{1}q_{1} \text e^{\mu _{1}(t-T)}+c_{2}q_{2} \text e^{\mu _{2}(t-T)}\right ] +{\bar{\Lambda }}_{2}, \end{aligned}$$

where \(\mu _{1}\) and \(\mu _{2}\) are the two roots of the quadratic equation

$$\begin{aligned} \mu ^{2}-\left( \alpha _{1}+\alpha _{2}\right) \mu -\alpha _{1}\alpha _{2}=0 \end{aligned}.$$

\(\left( 1,q_{1}\right) ^{\text T}\) and \(\left( 1,q_{2}\right) ^{\text T}\) are the appropriate eigenvectors associated with the eigenvalues \(\mu _{1}\) and \(\mu _{2}\) of the coefficient matrix of the adjoint ODE, and

$$\begin{aligned} {\bar{\Lambda }}_{1}=\frac{u_{\max }}{u_{\max }-\alpha _{1}},\ \ \ \ \ \bar{ \Lambda }_{2}=2{\bar{\Lambda }}_{1}. \end{aligned}$$

The two constants \(c_{1}\) and \(c_{2}\) are chosen so that adjoint end conditions \(\Lambda _{1}(T)=\Lambda _{1}(T)=0\). The switch point \(t_{\text s}\) is then determined by

$$\begin{aligned} \Lambda _{1}(t_{\text s})=\left [ c_{1} \text e^{\mu _{1}(t_{\text s}-T)}+c_{2} \text e^{\mu _{2}(t_{\text s}-T)}\right ] +\frac{u_{\max }}{u_{\max }+\alpha _{1}}=1. \end{aligned}$$

The switch condition \(\Lambda _{1}(t_{\text s})=1\) for determining the switch point \(t_{\text s}\) is a nonlinear equation for \(t_{\text s}\) and can be solved (iteratively) by any one of the equation solver available on math software such as Mathematica. Instead, a single pass algorithm will be formulated in the next subsection that would lead to an exact explicit solution for the problem and hence intrinsically more efficient (even if it is not necessarily a huge gain in efficiency for the present relatively simple problem).

6.2.3 A Single Pass Algorithm for the Switch Point

For a more efficient and effective algorithm for computing the switch point \(t_{\text s}\), we interchange the roles of t and \(\lambda _{1}\) and work with the latter as the independent variable. In that case, the two adjoint equations and the terminal conditions take the form

$$ \begin{aligned} \frac{\mathrm{{d}}t}{\mathrm{{d}}\lambda _{1}} & = \left( \frac{\mathrm{{d}}\lambda _{1}}{\mathrm{{d}}t}\right) ^{-1}= \frac{1}{\left( \alpha _{1}+u_{\max }\right) \lambda _{1}-\alpha _{1}\lambda _{2}-u_{\max }},\ \ \ \ \ \ t(0)=T, \\ \frac{\mathrm{{d}}\lambda _{2}}{\mathrm{{d}}\lambda _{1}} & = \frac{ {\frac{\mathrm{{d}}\lambda _{2}}{\mathrm{{d}}t}}} {\frac{{\mathrm{{d}}\lambda _{1}}}{\mathrm{{d}}t}}=\frac{\alpha _{2}\left( \lambda _{2}-2\lambda _{1}\right) }{\left( \alpha _{1}+u_{\max }\right) \lambda _{1}-\alpha _{1}\lambda _{2}-u_{\max }} ,\ \ \ \ \ \ \ \lambda _{2}(0)=0. \end{aligned} $$

This new system is de-coupled so that we may solve the second for \(\lambda _{2}=f(\lambda _{1})\) and then the first to get \(t=g(\lambda _{1})\). The switch point \(t_{\text s}\) is then given immediately by the end point value

$$\begin{aligned} t_{\text s}=g(1). \end{aligned}$$

Since the second equation for \(\lambda _{2}=f(\lambda _{1})\) is scale invariant, we may transform it into a single first-order ODE which is separable which yields an exact explicit solution \(\lambda _{2}=f(\lambda _{1})\). The first ODE becomes

$$\begin{aligned} t=T-\int _{0}^{\lambda _{1}}\frac{\mathrm{{d}}\lambda _{1}}{u_{\max }+\alpha _{1}f(\lambda _{1})-\left( \alpha _{1}+u_{\max }\right) \lambda _{1}}\equiv g(\lambda _{1}). \end{aligned}$$

The integral on the right is an increasing function of \(\lambda _{1}\) so that

$$\begin{aligned} t_{\text s}=g(1)<T \end{aligned}.$$

The single pass solution above has been applied to obtain \(t_{\text s}\) for different combinations of \(\left\{ u_{\max },\alpha _{k},T\right\}\) to investigate the possibility of growth of \(C(t)\equiv R(t)+D(t)\) after the onset of RB-to-EB conversion. For that we need to work with the equations of state (36) and the corresponding initial conditions characterizing the growth dynamics of the RB and DB populations.

6.3 RB Population Growth in the Conversion Phase

Given that all population change rates are linear in the unknowns, the optimal control is expected (and can be shown) to be also the bang-bang control (7) for \(u_{\max }\) sufficiently large [21]. Prior to the onset of RB-to-EB conversion at \(t_{\text s}\), we have \(u_\mathrm{{op}}(t)=0\) with the first two ODE in (36) summed to give

$$\begin{aligned} C^{\prime }(t)=\alpha _{2}D(t),\ \ \ \ \ C(0)=N \end{aligned}$$
(42)

for \(C(t)=R(t)+D(t)\) and \(R(0)+D\left( 0\right) =N\). Since \(D(t)>0\) for \(t>0\) , we have \(C^{\prime }(t)>0\) so that the total RB population increases with time at least up to the onset of the RB-to-EB conversion at time \(t_{\text s}\).

For \(t>t_{\text s}\), we have from (7) \(u_\mathrm{{op}}(t)=u_{\max }\) so that

$$\begin{aligned} C^{\prime }=-u_{\max }R+\alpha _{2}D,\ \ \ \ \ C(t_{\text s})=C_{\text s}, \ \ \ \ \ t_{\text s}<t\leqslant T \end{aligned}$$
(43)

with \(C_{\text s}=R_{\text s}+D_{\text s}\) where \(R_{\text s}\) and \(D_{\text s}\) are, by continuity of the population sizes, the values of R and D at the switch point \(t_{\text s}\) known from their solutions for \(u=0\) in \([0,t_{ \text s}]\). Immediately after the onset of conversion, we have

$$\begin{aligned} C^{\prime }\left( t_{\text s}{\small +}\right) =-u_{\max }R(t_{\text s})+\alpha _{2}D\left( t_{ \text s}\right) . \end{aligned}$$

We have then the following important conclusion pertaining to the purpose of our investigation.

Proposition 6

The total C(t) of RB and DB populations continues to increase at least for a period after \(t_{\text s}\) if and only if

$$\begin{aligned} u_{\max }R(t_{\text s})<\alpha _{2}D\left( t_{\text s}\right) . \end{aligned}$$
(44)

This proposition establishes the feasibility of the four-form model for the continual growth of the total reticulate bodies (RBs and DBs) after the onset of conversion within a range of the conversion capacity as stipulated by (44). It should be noted however that to close the model with only two types of reticulate bodies, R(t) and D(t), the characterization of the RB-to-EB conversion by (37) is rather artificial since only a fraction of the R(t) in the four-form model is eligible to convert. Furthermore, that portion changes with time as more and more RBs become size eligible. A mechanism will have to be developed to capture this feature before we should proceed with additional work on this model. In this context, the bound (44) on \(u_{\max }\) for continual growth of C(t) beyond \(t_{\text s}\) is likely to be overestimated, since only a fraction \(R^{(c)}(t)\) of R(t) is eligible to convert to EB (via IB or not) at any time t. The actual RB population after the onset of conversion is likely to be larger than that determined by our model, assuming no other factors affect the estimate.

This deficiency notwithstanding, the condition \(\alpha _{2}D_{\text s}>u_{\max }R_{\text s}\) of the proposition appears to be met by the data reported in [11]. Crude estimates using the data there lead to the comparisons in Table 3 showing \(\alpha _{2}D>u_{\max }R\) at least for the entire interval \(\left( 24,28\right)\) hpi with \(R(t)+D(t)\) begins to decline sometime before 32 hpi. With Proposition 6 pretty much met by the data on the infected cells observed [11], and the infected host cell of Fig. 4 qualitatively similar to the data reported on total RB sum C(t) of RB and DB, the four-form model of this section appears also capable of rectifying the deficiencies on the RB growth in the conversion phase of the simpler models of [8, 21] even if the model closure seems not biologically realistic.

Table 3 Growth rate estimates
Full size table

7 Summary of Findings and Implications

The comprehensive examination of chlamydial inclusions by three-dimensional electron microscopy reported in [11] provided a rich trove of data on the unusual Chlamydia developmental cycle. Among the developmental features of the bacterium’s reticulate bodies observed relevant to our investigation are

  1. (i)

    the alternative bacterial fate choices between division into two RBs, or conversion into an EB;

  2. (ii)

    the existence of a period of no conversion (a conversion holiday) in the first half of the developmental cycle;

  3. (iii)

    a putative small size threshold for the RB-to-EB conversion; and

  4. (iv)

    the RB size reduction through repeated divisions to reach this small conversion threshold size.

These unusual features prompted the initiation of theoretical research projects to gain some insight into the bacteria’s bio-theoretics of development through modeling, analysis, and computation. Among these projects are the formulation and analysis of optimal control models in [8, 21]. They enabled us to uncover the biological force that drives the developmental cycle by providing a link between the progression of the RB-to-EB conversion, including the conversion holiday at the start of the infection, observed empirically in [11], and the maximization of the spread of infection when the host cell lyses. In all cases, these models require the strategy of maximizing the terminal EB population (that are solely responsible for new infection) to be a bang-bang conversion rate, a requirement consistent with a conversion holiday observed at the start of all C trachomatis life cycles.

The optimal control models in [8, 21] are not sufficiently fine-grained for their predictions to reflect all actual features of the bacteria’s developments. While it is rather remarkable that these models succeed in replicating the unusual conversion holiday at the start of the life cycle, it is deficient in characterizing the evolution of the RB population after the onset of the RB-to-EB conversion at the switch point \(t_{\text s}\). The total (lumped) RB units predicted to decrease immediately after the onset of conversion are at variance with the data reported in [11]. As summarized in Table 1, the RB total (averaged over the inclusions observed) should increase for \(t>t_{\text s}\) at least for a substantial period comparable to the duration of the conversion holiday before it begins to decline. By incorporating the small size threshold eligibility to convert, we formulate herein an n-stage RB model in which RB units reduce in size after each division by binary fission as observed in [11] and are eligible to convert when they reach a threshold size after a number of divisions (estimated to be less than six). We show that the solution of the corresponding optimal control problem continues to require an optimal bang-bang control (and hence a conversion holiday at the start of the infection). In addition, the total of all RB subpopulations may now continue to increase beyond the switch point for a period before declining as shown in the data reported in [11]. An important feature of the optimal solution involves a bound on the maximum conversion capacity \(u_{\max }\) (in relation to other bacterial and host cell characteristics such as the natural growth rate constants of the various RB subpopulations of different vintages, the host cell lysis time, etc.) for the RB total to grow beyond the switch point \(t_{\text s}\). The same total would decline immediately after the onset of conversion if \(u_{\max }\) should exceed the limit imposed by the bound. An even more interesting result is the permissible magnitude of \(u_{\max }\) for continual growth beyond \(t_{\text s}\) increases with n, and the number of divisions of the RB before the conversion eligibility is met as shown by the graphs of the RB total in Figs. 2, 3, 4, and 5.

To the extent that a reasonable range of \(u_{\max }-\alpha _{n}\) is required to attain the temporal profile of the total RB growth in the infected host cells (at least of the type observed and reported in [11]), data on the RB units of three or more vintages (\(n>2\)) would be needed for further quantification of the n-stage model. Until data for higher RB vintages (than DB) become available, we take gratification for the n-stage model’s capability to grow the total RB beyond the switch point while correctly reproduce the optimal strategy for species survival.