Immunophysical analysis of corneal neovascularization: mechanistic insights and implications for pharmacotherapy

The cornea lacks adaptive immune cells and vasculature under healthy conditions, but is populated by both cell types under pathologic conditions and after transplantation. Here we propose an immunophysical approach to describe postoperative neovascularization in corneal grafts. We develop a simple dynamic model that captures not only the well-established interactions between innate immunity and vascular dynamics but also incorporates the contributions of adaptive immunity to vascular growth. We study how these interactions determine dynamic changes and steady states of the system as well as the clinical outcome, i.e. graft survival. The model allows us to systematically explore the impact of pharmacological inhibitors of vascular growth on the function and survival of transplanted corneas and search for the optimal time to initiatetherapy. Predictions from our models will help ongoing efforts to design therapeutic approaches to modulate alloimmunity and suppress allograft rejection.


TERMINOLOGY
We used the following convention for naming coefficients: for direct stimulation or suppression effect, we use α X1X2 in which X 1 is the actor parameter and X 2 is the target parameter. For example, α V B T is the contribution of blood vessels to trauma healing. If more than one parameter contributes to an action on another parameter, the corresponding coefficient contains all those parameters as actor parameters. Furthermore, for spontaneous decay, we used λ X for corresponding decay rate. Table II lists the exact definitions for all coefficients.

COEFFICIENT ESTIMATION
In order to estimate coefficients, we used a two-step process. First, we estimated the value of coefficients based on qualitative arguments and by generating results that agree closely with the experimental data. Then we used a quantitative approach to further improve our estimations.

Qualitative Estimation
To estimate each coefficient, we need to understand its contribution to the dynamics of the model and to compare it with the data from the dynamics of the real system [1]. Some coefficients were directly estimated and some others were obtained through screening reasonable ranges as follows: Figure 1B in [1] shows that value of V B for HR setting approaches 8; thus, we set V Bmax = 8. Besides, Figure 1C in [1] shows saturation of E in the case of HR at long time scales. Accordingly, the upper limit for E is close to 6. Simply by comparing experimental measurements for V B and V L (Figures 2 and 3 [1]) we can roughly estimate V Lmax = V Bmax /2 = 4. Adaptive immune response comes with a delay with respect to other components of our model. Consequently, we can ignore terms related to adaptive immunity for early times. As such, T only depends on V B (for t < 14 days). Figure 1B in [1] shows that for the LR case, the value of V B starts to fall at t = 7 days. This means that the value of T has changed significantly and reaches a low level at this time. Value of V B (t) 1+V B (t) during the first 7 days is roughly equal to 0.4 (based on simply averaging on this function for t = 0 and t = 7 and available data for V B ). We suppose that the value of T reaches 0.2 (recall T (0) = 1). By inserting all estimated values to equation 1 the modified version, dT (t) is easy to solve for the boundary conditions (T (0) = 1,T (7) = 0.2) and α V B T ≈ 0.5. On the other hand, we know that V B has a maximum at t = 7 days; thus, we have dV BN (t) dt | t=7 = 0. This relation leads to λ V B = α T V BN /35. For 0 < t < 7 we can assume T avarage = 0.6 and again V Baverage = 1.8. The modified version of equation 2, would be As a result, we have α T V BN = 1, λ V B = 0.03. V L in the LR case and at the early times has similar dynamics (equation 3). Value of V L in this regime is about half of the value of V B [1]( Figures 2B and 3B). Therefore, we can argue for α T V L = α T V BN /2. Besides, we know that the value of V B declines for LR. For large times, spontaneous decay is the only weakening term in the evolution of V B , then λ V L is the only responsible coefficient for this behavior and from results for long time and we found Average value of E in this range is equal to 0.85. Figure 1C in [1] shows that the behavior of E is analogous to the behavior of V B with half value. Here we use our estimations in previous parts to estimate values of the remaining coefficients. This After successfully creating the experimental results for early times, we consider longer time periods (t > 14 days). The rise of parameters after 14 days implies that the adaptive immune response has come to play prior to this time point. To see such behavior we have to set t d < 14 and it turns out that t d = 7 days fits best with experimental results. At long time scales, we know that V L reaches zero in the case of LR setting. As a result, λ E is the only coefficient with weakening effect on the dynamics of E and we found its value by monitoring the behavior of E in long time periods.
An accurate value of D is not available, however, we can come up with an estimate. The steady dynamics of other parameters in the long term suggest a steady dynamics for D itself. Thus, we expect a rise after a time delay and rather a constant value for D then. Therefore, D P should be set to satisfy this situation. Similar logical arguments can be made for D P . Its value should start to rise after t = t d and have smooth behavior in long times. We set coefficients in order to generate the mentioned behavior for D P . This rise and subsequent smooth behavior of D P leads to an initial rise and a subsequent steady dynamics in long time periods for D. Finally, experimental results were reproduced by tuning relevant parameters, while keeping the steady dynamics of D and D P invariant (see Figure S3). Again, we stress that our procedure involved an initial estimation which was subsequently tuned to reach the best fits with the experimental results. Figure S2 shows the qualitative agreement between our estimation and available experimental results.

Quantitative Refinement
The above mentioned qualitative estimations of the coefficients were conducted to find a reasonable domain for each coefficient, and for a more accurate fit we did a quantitative approach as we describe below: we consider three main parameters from LR model, i.e., evolution of V B , V L and E as our baseline experimental data. V B and E have been clinically measured but there exist no equivalent experimental data for V L . By comparing Figure 2B and Figure  3B in [1] we learn that V B V L > 2 for t < 28 days (for HR this ratio is different). Besides, for larger times, the value of V L decays ( Figure S4). As a result, we consider Table I data as our reference for fitting.
To obtain the best fit, we use the least square method and calculate R = i 56 t=0 (X i (t) sim − X i (t) exp ) 2 with X i = V B , V L and E for LR. Extreme minimum value of R as the best fit is obtainable by screening the space of coefficients and it turns out to be equal to 2.9 for the coefficient values listed in Table II. We note that these values  properly reproduce the behavior of HR setting as well.

SENSITIVITY ANALYSIS AND CONFIDENCE REGION
Here we perform sensitivity analysis to check the dependency of the model on our estimated parameter values. In other words, we ask how much deviation from each estimated value is needed to generate significant alteration in the dynamics of the model. We define a confidence region [2] around each coefficient value where the original dynamics of the model is conserved despite the introduced perturbation. To define the region, we consider R, as defined below, as the identification factor: We tune each coefficient until the value of R exceeds 4. Through this procedure we do a sensitivity analysis and obtain a confidence region for each coefficient (see Table II).

SHORTEST DURATION OF THERAPY
Since drugs typically carry some side effects, it is of high clinical value to find the shortest duration of therapy to minimize the side effects. To find the duration of therapy, we measured the time between the start of drug administration and the time when the edema score drops to below two. For both specific and non-specific approaches, we calculated the durations of therapy (see Figures S5 and S6). Note that contrary to Figure 9 of the article, here we do not ensure that after cessation of therapy edema score remains below two and therapy has not been administrated after 60 days.   VLmax The maximum value for lymphatic vessels (capacity of the tissue for lymphatic vessels) 4.6 2.

VBmax
The maximum value for blood vessels (capacity of the tissue for blood vessels) 8 0.8

Emax
The maximum value for edema (capacity of the tissue for edema) 5.35 2.2