Efficiency in Lung Transplant Allocation Strategies

TL;DR

This paper develops a new framework for lung transplant allocation that models patient and organ dynamics to optimize long-term patient survival, outperforming the current LAS system.

Contribution

It introduces a comprehensive stochastic model and an optimal allocation strategy that significantly improves patient life expectancy over existing methods.

Findings

Proposed strategy increases average total life by 7.7% compared to LAS.

Model accurately predicts optimal transplant timing based on patient health and organ availability.

Simulation results validate the effectiveness of the new allocation approach.

Abstract

Currently in the United States, lung transplantations are allocated to candidates according to the candidates' Lung Allocation Score (LAS). The LAS is an ad-hoc ranking system for patients' priorities of transplantation. The goal of this study is to develop a framework for improving patients' life expectancy over the LAS based on a comprehensive modeling of the lung transplantation waiting list. Patients and organs are modeled as arriving according to Poisson processes, patients' health status evolving a waiting time inhomogeneous Markov process until death or transplantation, with organ recipient's expected post-transplant residual life depending on waiting time and health status at transplantation. Under allocation rules satisfying minimal fairness requirements, the long-term average expected life converges, and its limit is a natural standard for comparing allocation strategies. Via the Hamilton-Jacobi-Bellman equations, upper bounds for the limiting average expected life are derived as a function of organ availability. Corresponding to each upper bound is an allocable set of (time, state) pairs at which patients would be optimally transplanted. The allocable set expands monotonically as organ availability increases, which motivates the development of an allocation strategy that leads to long-term expected life close to the upper bound. Simulation studies are conducted with model parameters estimated from national lung transplantation data. Results suggest that compared to the LAS, the proposed allocation strategy could provide a 7.7% increase in average total life. We further extended the results to the the allocation and matching of multiple organ types.

Figures

Tables

Table 1. Table 1: Means and standard deviations (SD) of averaged life outcomes (in days) with different allocation strategies over 200 200 200 independent simulations.

	Average Total Life (SD)	Average Life in Waiting List (SD)	Average Life Post-Transplant (SD)
Proposed	1839 (19)	489 (21)	1349 (20)
LAS	1708 (20)	554 (20)	1155 (15)
Refined LAS	1723 (21)	561 (21)	1162 (16)
Randomized	1612 (26)	469 (21)	1143 (30)

Table 2. Table 2: Comparison of Average Waiting Time with Different Allocation Strategies

	Average waiting time (un-transplanted) (SD)	Average waiting time (transplanted) (SD)
Proposed	234 (14)	256 (18)
LAS	267 (14)	286 (17)
Refined LAS	242 (15)	318 (20)
Randomized	213 (17)	257 (22)

Table 3. Table 3: Coefficient Estimates for In-Waiting-List Hazard for Death

Covariate	Group	Coefficient	$p$-value111Significance level: *** $p$-value $\le 0.001$, ** $p$-value $\le 0.01$, * $p$-value $\le 0.05$, . $p$-value $\le 0.1$.
Age (year)	3	0.014	(***)
Body mass index (BMI; $kg/m^2$)	1	0.116 $\times$ (20 - BMI) for BMI less than 20 $kg/m^2$	(***)
Ventilation status if candidate is hospitalized	1	-0.444 if continuous mechanical ventilation needed	$6.7\times {10}^{-14}$ (***)
Creatinine (serum, mg/dL)	1	0.220 if at least 18 years of age	(***)
Diabetes	1	0.173	(***)
Diagnosis Group A	3	0	NA
Diagnosis Group B	3	0.794	(***)
Diagnosis Group C	3	1.126	(***)
Diagnosis Group D	3	0.163	$3.3\times {10}^{-4}$ (***)
Detailed diagnosis: Bronchiectasis (Diagnosis Group A only)	3	0.182	$1.4\times {10}^{-3}$ (***)
Eisenmenger’s syndrome (Diagnosis Group B only)	3	-1.04	(***)
Lymphangioleiomyomatosis (Diagnosis Group A only)	3	-0.961	$1.7\times {10}^{-13}$ (***)
Obliterative bronchiolitis (Diagnosis Group D only)	3	-0.416	$8.3\times {10}^{-5}$ (***)
Pulmonary fibrosis, not idiopathic (Diagnosis Group D only)	3	0.014	0.70
Sarcoidosis with PA mean pressure greater than 30 mm Hg (Diagnosis Group D only)	3	-0.44	(***)
Sarcoidosis with PA mean pressure of 30 mm Hg or less (Diagnosis Group A only)	3	0.613	(***)
Forced vital capacity (FVC)	2	0.188 $\times$ (80 - FVC)/10 if FVC is less than 80% for Diagnosis Group D	(***)
Functional Status	1	-0.287 if no assistance needed with activities of daily living	(***)
Oxygen needed to maintain adequate oxygen saturation (80% or greater) at rest (L/min)	2	0.111 for Group B, 0.108 for Groups A, C, and D	(***)
PCO2 (mm Hg): current	1	0.222 if PCO2 is at least 40 mm Hg	(***)
PCO2 increase of at least 15%	1	-0.232 if PCO2 increase is at least 15%	$9.4\times {10}^{-15}$ (***)
Pulmonary artery (PA) systolic pressure (10 mm Hg) at rest, prior to any exercise	2	0.003 for Group A if the PA systolic pressure is greater than 40 mm Hg, 0.016 for Groups B, C, and D	(***), $9.3\times {10}^{-3}$ (***)
Six minute walk distance (feet) obtained while the candidate is receiving supplemental oxygen required to maintain an oxygen saturation of 88% or greater at rest. Increase in supplemental oxygen during this test is at the discretion of the center performing the test.	1	-0.075 $\times$ Six-minute-walk distance/100	(***)

Table 4. Table 4: Coefficient Estimates for Post-Transplantation Hazard for Death

Covariate	Group	Coefficient	$p$-value222Significance level: *** $p$-value $\le 0.001$, ** $p$-value $\le 0.01$, * $p$-value $\le 0.05$, . $p$-value $\le 0.1$.
Age (years)	4	$4.4\times {10}^{-3}\times$ (age - 45) if greater than 45 years of age	$2.1\times {10}^{-5}$ (***)
Creatinine (serum) at transplant (mg/dL)	1	0.177 if candidate is at least 18 years old	0.003 (**)
Creatinine increase of at least 150%	1	0.570 if increase in creatinine is at least 150%, and the higher value determining this increase is at least 1 mg/dL	0.422
Ventilation status if candidate is hospitalized	1	-0.05 if continuous mechanical ventilation needed	0.710
Diagnosis Group A	4	0	NA
Diagnosis Group B	4	0.263	0.017 (*)
Diagnosis Group C	4	0.268	0.004 (**)
Diagnosis Group D	4	0.171	0.009 (**)
Detailed diagnosis: Bronchiectasis (Diagnosis Group A only)	4	0.191	0.168
Eisenmenger’s syndrome (Diagnosis Group B only)	4	0.745	0.297
Lymphangioleiomyomatosis (Diagnosis Group A only)	4	-0.625	0.052 (.)
Obliterative bronchiolitis (not-retransplant, Diagnosis Group D only)	4	0.035	0.866
Pulmonary fibrosis, not idiopathic (Diagnosis Group D only)	4	-0.150	0.077 (.)
Sarcoidosis with PA mean pressure greater than 30 mm Hg (Diagnosis Group D only)	4	-0.230	0.078 (.)
Sarcoidosis with PA mean pressure of 30 mm Hg or less (Diagnosis Group A only)	4	-0.043	0.801
Oxygen needed to maintain adequate oxygen saturation (80% or greater) at rest (L/min)	4	$6.6\times {10}^{-3}$ for Group A $1.1\times {10}^{-3}$ for Groups B, C, and D	0.590, 0.845
Functional Status	1	-0.206 if no assistance needed with activities of daily living	0.001 (***)
Six minute walk distance (feet) obtained while the candidate is receiving supplemental oxygen required to maintain an oxygen saturation of 88% or greater at rest. Increase in supplemental oxygen during this test is at the discretion of the center performing the test.	1	$3.0\times {10}^{-4}\times$ (1200-Six minute walk distance), 0 if six-minute-distance-walked is at least 1200 feet	$6.2\times {10}^{-10}$ (***)

Equations

$${\tilde{X}}^{(k)}(s)=X^{(k)}(s)\cdot I(s<T_T^{(k)})+(n+1)\cdot I(s\ge T_T^{(k)}),$$

$$\mathcal{F}_{t}=\sigma\big{(}T_{k},S_{i},a_{i},\{X^{(k)}(s):s\in[0,t-T_{k}]\}:T_{k},S_{i}\in[0,t]\big{)}.$$

$$W_t=\{(t-T_k,{\tilde{X}}^{(k)}(t-T_k)):T_k\le t,{\tilde{X}}^{(k)}(t-T_k)\notin \{0,n+1\}\}.$$

$$\mathcal{G}_{t}=\sigma\big{(}T_{k},S_{i},\{X^{(k)}(s):s\in[0,t-T_{k}]\}:T_{k},S_{i}\in[0,t]\big{)}\vee\sigma(a_{i}:S_{i}<t),$$

$$W_t^{\prime }=\{(t-T_k,X^{(k)}(t-T_k)):T_k\le t,k\notin \{a_i:S_i<t\},X^{(k)}(t-T_k)\ne 0\}.$$

$$P((S_j-T_{a_j},X^{(a_j)}(S_j-T_{a_j}))=x\mid {\mathcal{G}}_{S_j})=\Gamma (x,W_{S_j}^{\prime }).$$

$$\begin{gathered} P((S_j-T_{a_j},X^{(a_j)}(S_j-T_{a_j}))=x \\ \mid W_{S_j}^{\prime },\{S_i,T_k:i,k\in \mathbb{N}\},\{X^{(k)}(s),R^{(k)}(s):s\ge 0,k\in \mathbb{N}\}) \\ =\Gamma (x,W_{S_j}^{\prime }). \end{gathered}$$

$$\underset{t\ge 0}{\inf }\{X^{(k)}(t)=0\}\wedge T_T^{(k)}+\sum _{i\in \mathcal{X}}{\int }_0^T{R}_i^{(k)}(s)d{\tilde{N}}_{i,n+1}^{(k)}(s),$$

$$\sum _{i=1}^n{\int }_0^T{\tilde{R}}^{(k)}(s)d{\tilde{N}}_{i,n+1}^{(k)}(s)+\underset{t\ge 0}{\inf }\{X^{(k)}(t)=0\},$$

$$\lim_{t\rightarrow\infty}\frac{1}{N_{t}}E\Bigg{[}\sum_{j=1}^{O_{t}}\tilde{R}^{(a_{j})}(S_{j}-T_{a_{j}})\Bigg{]},$$

$${\pi }_i(s)=\underset{t\to \infty }{\lim }\frac{1}{N_t}\sum _{k=1}^{N_t}I({\tilde{X}}^{(k)}(s)=i).$$

$$d{\mathbf{\Psi }}_s={\text{uppi}}_{s-}{\mathbf{Q}}_{s}ds-d{\text{uppi}}_s,$$

$$\lim_{t\rightarrow\infty}\frac{1}{N_{t}}{E}\Bigg{[}\sum_{j=1}^{O_{t}}\tilde{R}^{(a_{j})}(S_{j}-T_{a_{j}})\Bigg{]}=\int_{0}^{T}d\boldsymbol{\Psi}_{s}\tilde{\boldsymbol{\upmu}}_{s}.$$

$${\int }_{[0,T]}d{\mathbf{\Psi }}_s{1}_n\le \frac{\rho }{\tau },$$

$${\int }_{[0,T]}d{\mathbf{\Psi }}_s{\widetilde{\text{upmu}}}_s-c{\int }_{[0,T]}d{\mathbf{\Psi }}_s{1}_n$$

$${\Psi }_i^c(\{s\})=I({\phi }_i^c(s)>c)\cdot {\pi }_i(s-)$$

$$d{\Psi }_i^c(s)=I({\phi }_i^c(s)>c)\cdot \sum _{j\in \mathcal{X}}{\pi }_j(s-)q_{ji}(s)ds$$

$${\phi }_i^c(s)=\frac{{\tilde{\mu }}_i(s)-{\eta }_i^c(s+)}{1-{\gamma }_i^c(s+)},$$

$${\text{upeta}}^c(s+)={\int }_{(s,T]}{\text{Prodi}}_{(s,t)}(\mathbf{I}+d\mathbf{Q})(\mathbf{I}-d{\mathbf{\Lambda }}^c)\cdot d{\mathbf{\Lambda }}_t^c{\widetilde{\text{upmu}}}_t,$$

$${\text{upgamma}}^c(s+)={\int }_{(s,T]}{\text{Prodi}}_{(s,t)}(\mathbf{I}+d\mathbf{Q})(\mathbf{I}-d{\mathbf{\Lambda }}^c)\cdot d{\mathbf{\Lambda }}_t^c{1}_n,$$

$${\phi }_i^c(s)=\frac{{\tilde{\mu }}_i(s)-{\eta }_i^c(s+)}{1-{\gamma }_i^c(s+)}$$

$${\int }_{[0,T]}d{\mathbf{\Psi }}_s^{c_1}{1}_n\le {\int }_{[0,T]}d{\mathbf{\Psi }}_s^{c_2}{1}_n.$$

$${\int }_{[0,T]}d{\mathbf{\Psi }}_s^{c_1}{\widetilde{\text{upmu}}}_s\le {\int }_{[0,T]}d{\mathbf{\Psi }}_s^{c_2}{\widetilde{\text{upmu}}}_s.$$

$${\mathcal{A}}_c=\{(i,s)\in \mathcal{X}\times [0,T]:{\phi }_i^c(s)>c\},$$

$${\mathcal{A}}_{c_1}\subseteq {\mathcal{A}}_{c_2}.$$

$$\sum _{w=1}^{N_o}{\int }_0^{T}d{\mathbf{\Psi }}_s^{(w)}{\widetilde{\text{upmu}}}_s^{(w)}.$$

$${\int }_{[0,T]}d{\mathbf{\Psi }}_s^{(w)}{1}_n\le \frac{{\rho }^{(w)}}{\tau },$$

$$\sum _{w=1}^{N_o}{\Psi }_i^{(w)}(\{s\})\in [0,{\pi }_i(s-)]\text{ and }\sum _{w=1}^{N_o}d{\Psi }_i^{(w)}(s)\in [0,\infty ),\text{ if }{\pi }_i(s-)\ne 0,$$

$$\sum _{w=1}^{N_o}d{\Psi }_i^{(w)}(s)\in [0,\{{\mathbf{\pi }}_{s-}d{\mathbf{Q}}_s{\}}_i],\text{ if }{\pi }_i(s-)=0,$$

$$\sum _{w=1}^{N_o}\left({\int }_{[0,T]}d{\mathbf{\Psi }}_s^{(w)}{\widetilde{\text{upmu}}}_s-c^{(w)}{\int }_{[0,T]}d{\mathbf{\Psi }}_s^{(w)}{1}_n\right),$$

Full text

Lung transplantation,

Lung transplantation score (LAS),

Allocation rules,

Markov decision models,

Hamilton-Jacobi-Bellman equations,

Simulation studies,

,

and

1 Introduction

1.1 The Lung Transplantation Score (LAS)

Lung transplantation was first performed in 1963 and has evolved over the years to become a standard treatment for patients with advanced lung diseases (Kotloff and Thabut (2012)). There are not enough organs available for transplantation to meet the need of patients waiting for lungs, and efficient allocation of lung transplantations to patients is vital.

Currently in the United States, patients in need of lung transplantation are registered to waiting lists managed by the United Network for Organ Sharing (UNOS). As lungs become available for transplant, they are allocated to candidates in waiting lists on the basis of age, geography, blood type (ABO) compatibility, and the Lung Allocation Score (LAS). The LAS was first implemented in 2005, aiming to reduce the risk of in-waiting-list mortality and to prolong post-transplant residual life of organ recipients (Valapour et al. (2017)). The LAS is calculated as $100\cdot(\textrm{PTAUC}-2\cdot\textrm{WLAUC}+730)/1095$, where WLAUC is the estimated in-waiting-list life expectancy during an additional year if not transplanted and PTAUC is the estimated post-transplant life expectancy during the first year. Both in-waiting-list and post-transplant survivals are computed using proportional hazards models given current values of relevant characteristics of patients. The expected numbers of days of survival are calculated by integrating the areas under the covariate-specific survival curves within the first year, under the assumption that a patient’s health status would remain constant. When organs become available, patients in the waiting list are ranked according to their LAS values. Candidates with higher LAS values are given higher priorities for transplantation (UNOS (2017)).

The LAS is recognized to be imperfect. In February 2015, the Organ Procurement and Transplantation Network (OPTN) managed by UNOS implemented a revision of the LAS. Modifications were made to the covariates and their weights in the calculation to better predict patients’ survival (Valapour et al. (2017)).

There are other concerns about the LAS besides the selection of covariates. Studies have demonstrated that the emphasis on reducing the risk of in-waiting-list mortality in the LAS calculation may have the effect of increasing post-transplant mortality, since patients with the highest LAS are often in the worst conditions and subject to high post-transplant risk of death (Liu et al. (2010), Russo et al. (2010, 2011), Merlo et al. (2009)). And hence the question of whether the patients prioritized by the LAS are those who would benefit the most from transplantation. Studies also suggested that the LAS focused on the 1-year survival and failed to assess the long-term benefit of lung transplantation for the recipients. Results in Maxwell et al. (2014) showed that the 5-year survival had in fact slightly decreased after the implementation of the LAS.

1.2 Previous Studies on Organ Transplantation

Organ allocation has been studied from the policy maker’s perspective. The goal is generally characterized as assignments of organs to patients in the waiting list that optimize expected outcomes. Early works can be traced back to Derman, Lieberman and Ross (1971, 1975), Albright and Derman (1972) and Righter (1989), in which the authors solved for the optimal solution for allocating sequentially arrived resources to a finite number of subjects. The focus of these studies was on general resource allocation problems and for this purpose, the modeling was simplified and not tailored specifically to the organ transplantation problem.

Recent studies have involved increasing model complexity in studying the organ allocation problem. Zenios (1999) and Zenios, Wein and Chertow (1999); Zenios, Chertow and Wein (2000) modeled the waiting list for kidney transplantation with a deterministic fluid model: patients and organ donors were categorized into classes based on their demographic, immunological, and physiological characteristics and different categories of patients and organs flow in and out of the waiting list at class-specific rates. An optimal allocation rule was derived under the standard of a linear combination of quality-adjusted life expectancy and equity. Akan et al. (2012) modeled the waiting list for donated livers with a fluid model, in which patients were categorized into multiple classes according to their health status. Patients in one class were assumed to be able to flow only into adjacent categories with deterministic rates. Optimal allocation strategies were developed under a standard combining expected total number of in-waiting-list deaths and quality-adjusted total life expectancy.

Organ allocation had also been studied from an individual patient’s perspective by, for example, Ahn and Hornberger (1996), Hornberger and Ahn (1997), Su and Zenios (2004, 2005, 2006), Alagoz et al. (2004, 2007a, 2007b)), and Ata, Ding and Zenios (2018). These studies focused on optimal patient strategies for accepting or rejecting offers of organs.

Aspects of the organ transplantation problem other than organ allocation have also been studied. For example, the impact of available cadaveric kidneys on the number of candidates in the waiting list (Ruth, Wyszewianski and Herline (1985)), kidney exchange programs (Roth, Sönmez and Ünver (2004); Roth, Sönmez and Utku Ünver (2005), Ashlagi et al. (2011, 2012), Cechlárová and Lacko (2012)), the design of policies that meet fairness constraints chosen by the policy maker (Bertsimas, Farias and Trichakis (2013)), and comparing the efficiency and fairness of three proposed heart allocation policies via simulation studies (Hasankhani and Khademi (2017)). In one of the few studies on lung transplantation, Vock et al. (2013) evaluated the survival benefit of lung transplantation with the LAS from a causal inference perspective.

1.3 Overview of Model and Results

Relative to the richness of the literature on kidney and liver transplantation, few attempts have been made to study organ allocation in lung transplantation. Yet models for kidney and liver transplantation waiting lists cannot be appropriated wholesale in the service of modeling lung transplantation as several characteristics distinguish the transplantation of lungs from the transplantation of other organs. Here we propose a model for the lung transplantation waiting list based on common characteristics shared by different organ types and on the characteristics unique to lungs.

A waiting list for lung transplantation consists of patients added at random calendar times with different characteristics such as diagnosis, demographics, and health indicators. Here patient arrivals are modeled by a calendar time homogeneous Poisson process. Patients’ initial states upon arrival are modeled as sampled randomly from a finite state space.

Patients’ states, especially their health status, may change during their sojourn on the waiting list. As patients’ states change, their hazards for in-waiting-list death and expected post-transplant residual life change accordingly, which in turn affect patients’ priorities for transplantation. Given the significant implications of changes in patients’ states on allocation decisions, it is crucial to model the trajectory of patients’ states accurately and comprehensively.

Here we model counterfactual patient health status trajectories, that is, trajectories would be observed without transplantations, as independent and identically distributed continuous-time Markov processes. The model proposed here reflects UNOS practice and that patients evolve in continuous time, and differs from those in the previous studies, in that a patient’s health status trajectory is indexed by the waiting time since the patient arrives, and the change in health status is assumed to be inhomogeneous with respect to waiting time. Waiting time since listing is an important indicator for a patient’s functional age, which is increasingly accepted in predicting potential outcomes (Kotloff and Thabut (2012)). The inhomogeneous assumption enables the model to capture changes in the health status transition rate as waiting time increases. A scenario that would benefit from the proposed modeling, for instance, would be when patients experience accelerated deterioration in health after waiting for a long period in a serious condition. In addition, instead of only allowing transitions between adjacent states, the proposed model allows the consideration of more complex patterns of transitions in patients’ health status.

Organs, either from deceased or living donors, also become available at random calendar times. Here organ arrivals are modeled as a homogeneous Poisson process that is independent of patient arrivals and transitions in health status. When an organ becomes available, a patient in the waiting list is selected for transplantation. Due to the short time between when an organ becomes available and when the organ is no longer viable for transplantation, available organs are transplanted almost instantly to the selected patient. An organ recipient’s health status at the time of transplantation affects the post-transplantation residual life (UNOS (2011)). Candidates who are not transplanted will remain in the waiting list for future transplantation opportunities. Except in rare circumstances, patients leave the waiting list only in case of transplantation or death.

In contrast to the case of kidney and liver transplantation, living donor transplantation is extremely rare in lung transplantation: according to UNOS (2011), 8,674 patients were in the waiting list from year 2009 to 2011, 5,172 received deceased donor transplant and only two received living donor transplant. Therefore our focus here is on deceased donor transplantation only. Statistics in UNOS (2011) also show that from 2009 to 2011 only 20 patients out of 5,192 refused the transplantation offer, indicating patient choice is practically negligible in modeling lung transplantation. Therefore, we ignore the possibility of organ refusal and focus on the policy maker’s perspective in studying lung allocation.

The policy maker uses allocation rules to determine which patient in the current waiting list is selected for transplantation when an organ is available. Allocation rules may make use of any information of the current state or history of the patients in the waiting list. There are ethical issues to be considered regarding fairness or equity. Zenios, Chertow and Wein (2000) defined absolute equity with which the discrepancies in the outcome across patient groups would be eliminated and relative equity in which the first-come-first-transplanted policy was used as the golden standard. Su and Zenios (2006) and Ata, Ding and Zenios (2018) defined equity as equal chances of receiving an organ for all patient groups. Akan et al. (2012) suggested to remove factors such as age, gender and race in making allocation decisions.

In our study, allocation rules that are not consistent over time or that are affected by factors other than waiting times and relevant health characteristics of patients are viewed as inequitable. Relevant health characteristics are those that predict outcomes by which allocation rules are compared, such as the current and future hazards for death and/or post-transplantation survival. Fair allocations should depend on waiting times and health characteristics of the patients in the waiting list at the time of organ arrivals and possible independent randomizations. Allocation rules obtained by defining an index as a function of a patient’s waiting time and relevant characteristics and allocating to the patient in the list with the highest index value (or randomly among patients tied for the highest value) are guaranteed to be fair.

Some standards for comparing allocation rules that have been considered may be characterized in terms of expectations. Defining a standard in terms of an expectation, however, does not unambiguously define an optimization problem. Such a definition begs the question of which patients’ expectations are to be maximized. Implicitly, such standards refer to the long-term averages of life quantities of patients in the waiting list as time goes to infinity.

Here we take the expected total life (in-waiting-list life plus post-transplant life) averaged over all patients as calendar time increases as the metric for evaluating allocation strategies. We show that under the minimal fairness constraints, the waiting list has a unique limiting distribution. On average, patients entering the waiting list are transplanted according to a unique allocation-rule-specific limiting transplantation rate, which is a function of waiting time and health status. The long-term average of patients’ total life exists and is a functional of the limiting transplantation rate. The method used here can be extended readily to include the long-term average of essentially any aspect of patient expectation.

It is also shown that the expected proportion of transplanted patients is bounded by the ratio of the intensity of organ arrivals to that of patient arrivals, and that the transplantation rate satisfies boundedness constraints related to the counterfactual transition rates. We begin the search for the optimal fair allocation rule by solving for the limiting waiting time and health status specific transplantation rate that optimizes the long-term average life subject to the boundedness constraints.

The Hamilton-Jacobi-Bellman equations are used to characterize the form of the optimal limiting transplantation rate. The optimal rate is nonzero whenever the difference between the expected residual life with an immediate transplant and the expected residual life without an immediate transplant is greater than one minus the probability of future transplantations, scaled by a penalty parameter associated with the ratio of organ arrival rate to patient arrival rate.

Not every transplantation rate satisfying the constraints is a limiting rate for some allocation strategy. In particular, the rate corresponding to the optimal solution may not be achievable, as patients in the waiting list that are alive and not transplanted do not necessarily include those at waiting times and states where the optimal rate is non-zero when an organ arrives.

Here we propose an allocation strategy that is designed so that the corresponding waiting time and health status specific transplantation rate is close to the rate given by the optimal solution. The proposed allocation strategy relies on a critical monotonicity property: as the penalty parameter associated with organ availability decreases, combinations of waiting time and health state become allocable with non-zero optimal transplantation rate in a monotonic manner, and the order of each combination becoming allocable makes a natural index of transplantation priority.

A comprehensive simulation study is conducted with model parameters estimated from the national lung transplantation data provided by UNOS to examine the effect of the proposed strategy comparing to the current LAS system. Results suggest that the proposed allocation strategy can provide a gain of at least 7% in expected average total life relative to the LAS.

We further extend the single-organ-type results to the case of multiple organ types. The optimal transplantation rates for all organ types are derived and an allocation strategy to match and allocate organs of each type to patients is developed based on monotonicity properties of the optimal allocation rates.

We would also like to note that although our focus in this article is on lung transplantation, the model and methods developed can be applied to a wide range of applications, including the transplantation and matching of other organ types and allocations of scarce (medical) resources to candidates in complex transitions of states in general.

In what follows, Section 2 introduces notation in modeling the lung transplantation waiting list. Section 3 characterizes the comparison of allocation rules in terms of a constrained optimization problem and states the optimal solution. Section 4 proposes a practical allocation strategy in terms of an allocation index based on the optimal solution. Section 5 extends the results to the multiple-organ-type case. Section 6 compares average total life of patients with the proposed allocation strategy and the LAS in simulation studies. Section 7 provides further numerical results of the multiple organ type scenario. Proofs of the theoretical results are postponed to the Appendix and Supplementary Material.

2 The Waiting List

Let $\tau<\infty$ denote the intensity of the Poisson process of patient arrivals to the waiting list, and let $0<T_{1}<T_{2}<\dots$ denote the patients’ arrival times. For convenience, let $N_{t}$ denote the number of patient arrivals up to calender time $t$. Let $\rho<\infty$ denote the intensity of the Poisson process of organ arrivals, and let $0<S_{1}<S_{2}<\dots$ denote the arrival times of organs. Let $O_{t}$ denote the number of organ arrivals up to calendar time $t$. Assume organ arrivals are independent of patient arrivals. We are interested in settings when $\rho<\tau$, as otherwise the supply of organs would generally meet the demand, and there would be no need for an allocation rule.

Let $\mathcal{X}=\{0,1,...,n\}$ denote the finite set of possible patient health status, in which [math] denotes the absorbing state corresponding to death. Let $\{X(s):s\geq 0\}$ denote a generic health status trajectory indexed by waiting time since arrival. Suppose $X(s)$ is sufficiently detailed that $\{X(s):s\geq 0\}$ is a càdlàg Markov process. Denote the transition kernel of the process by $p_{s,t}(i,j)$, and the infinitesimal generator by $q_{ij}(s)=\lim_{t\downarrow s}(p_{s,t}(i,j)-I{(i=j)})/(t-s)$, where $I{(\cdot)}$ denotes the indicator function. Relative to the granularity of measurements of patient health status, there are no patient states that lead inevitably to sudden transitions, so it may safely be assumed that the max total transition rate, $\sum_{j\neq i}q_{ij}(s)<\infty$. In addition the expected number of transitions in any finite waiting time interval is finite, for which a sufficient condition might be $\int_{0}^{T}q_{ij}(s)ds<\infty$ for any $i,j\in\mathcal{X}$. Let the vector $\boldsymbol{p}=(p_{1},\dots,p_{n})$ denote the distribution of patients’ initial states upon arrival, where $p_{i}=P(X(0)=i)$ for each $i\in\mathcal{X}$. As patients must enter the waiting list alive, $\sum_{i=1}^{n}p_{i}=1$. Let $T$ be the upper bound of a patient’s waiting time, i.e., $P(X(s)=0)=1$ for all $s\geq T$.

The effect of allocation to a patient is measured by the difference between the patient’s expected residual life with a transplant and expected residual life without a transplant. Let $R(s)$ denote a generic post-transplant residual life if transplantation occurs at waiting time $s\geq 0$. Assume that the characterization of patient health states is sufficiently informative such that for any $s\geq 0$, $R(s)$ and $\sigma(\{X(u):u<s\})$ are conditionally independent given $X(s)$. Since post-transplant residual life is non-negative and bounded, $R(s)\geq 0$ and $\sup_{s\geq 0,i\in\mathcal{X}}E\left(R(s)\mid X(s)=i\right)<\infty$.

Let $\{(X^{(k)}(s),R^{(k)}(s)):s\geq 0\}$ denote the counterfactual health status trajectories and post-transplant residual life processes of the $k$th arriving patient. $\{(X^{(k)}(s),R^{(k)}(s)):s\geq 0,\,k\in\mathbb{N}\}$ are independent and identically distributed copies of the generic pair $\{(X(s),R(s)):s\geq 0\}$, and are independent of $\{T_{i},S_{j}:i,j\in\mathbb{N}\}$.

An allocation sequence $\{a_{j}:j=1,2,\dots\}$ is a random sequence of patient indices, where $a_{j}$ is the index of the recipient of the $j$th organ. Organs are allocated promptly after being retrieved to keep their functionality, and the short delay in allocation after organ arrival is omitted. For each $j$, $T_{a_{j}}\leq S_{j}$, $X^{(a_{j})}(S_{j}-T_{a_{j}})\neq 0$, and $a_{j}\neq a_{k}$ if $j\neq k$, which reflects that patients who receive transplantation must have entered the waiting list before the organ and must be alive at the moment of allocation, and that patients exit the waiting list upon receiving a transplant.

Let $T_{T}^{(k)}$ denote the $k$th patient’s waiting time at transplantation, that is, $T_{T}^{(k)}=S_{j}-T_{k}$ if $a_{j}=k$ and $T_{T}^{(k)}=\infty$ if $k\notin\{a_{j}:j\in\mathbb{N}\}$. Let $(n+1)$ denote the post-transplant state and let $\mathcal{X}_{0}=\mathcal{X}\cup\{n+1\}$ denote the augmented patient state space. Denote the $k$th patient’s actual trajectory by $\{\tilde{X}^{(k)}(s):s\geq 0\}$, so that

[TABLE]

and $\tilde{X}^{(k)}(s)$ is a càdlàg process on $\mathcal{X}_{0}$.

To formulate the waiting list, we first define a filtration that describes the information relevant to events in the waiting list up to time $t$, including patient and organ arrivals, patients’ counterfactual health state transitions and allocation decisions:

[TABLE]

At any time, the current waiting list consists of waiting time and health status of patients who have arrived and have not died nor transplanted. Let $w_{0}$ denote the state of no patient in the list. On $\mathcal{W}=\cup_{d=1}^{\infty}([0,T]\times\mathcal{X})^{d}\cup\{w_{0}\}$ and its Borel $\sigma$-algebra $\mathcal{B}(\mathcal{W})$, define the waiting list process

[TABLE]

$\{W_{t}:t\geq 0\}$ is adapted to $\{\mathcal{F}_{t}\}$. Since discontinuities of $\{W_{t}:t\geq 0\}$ can only be a result of patient or organ arrivals or transitions in health status, which are all càdlàg, $W_{t}$ is càdlàg.

We define another filtration generated by events in the waiting list up to time $t$, but excluding events associated with the allocation decision at $t$.

[TABLE]

Denote the waiting list before potential allocation at calendar time $t$ by

[TABLE]

For any $j\in\mathbb{N}$, $S_{j}$ is a $\{\mathcal{G}_{t}\}$-stopping time. Let $\mathcal{G}_{S_{j}}=\{A:A\cap\{S_{j}\leq t\}\in\mathcal{G}_{t},\,\forall\,t\geq 0\}$, then $W^{\prime}_{S_{j}}$ is measurable with respect to $\mathcal{G}_{S_{j}}$ as $W^{\prime}_{s}$ is progressively measurable.

Allocation sequences should satisfy fairness requirements. Our definition of fairness requires that given the waiting list up to the moment of allocation, the choice of the organ recipient should only depend on waiting time and health status of patients who are currently alive and not transplanted, and a randomization that is conditionally independent, given the current state of the waiting list, of the history of the waiting list. In formal terms,

Definition 2.1**.**

An allocation sequence $\{a_{1},a_{2},\dots\}$ is termed fair, if there exists a function $\Gamma(\cdot,\cdot)$, such that

for any $x\in[0,T]\times\mathcal{X}$, $\Gamma(x,\cdot)$ is a measurable function on $(\mathcal{W},\mathcal{B}(\mathcal{W}))$, and for any $w\in\mathcal{W}\backslash\{w_{0}\}$, $\Gamma(\cdot,w)$ a probability measure on $\{x:x\in w\}$, 2. 2.

for any $j\in\mathbb{N}$ and $A\in\mathcal{B}(\mathcal{W})$,

[TABLE]

In the above definition, $\Gamma$ is the allocation rule that determines the probabilities of allocating an available organ to patients in the current waiting list. Note the definition implies that given the waiting times and health status of patients in the current waiting list, the allocation probabilities are invariant to patient indices and consistent with respect to the calendar time.

Moreover, only realistic fair allocation rules that are non-informative of patients’ post-transplant residual life and future events in the waiting list are considered, and thus

[TABLE]

The effect of an allocation rule $\Gamma$ is measured by the long-term average expected total life, including life in waiting list and post-transplantation, of all patients ever enter the waiting list. The $k$th patient’s total life is equal to

[TABLE]

where $d\tilde{N}_{i,n+1}^{(k)}(s)=1$ if $T^{(k)}_{T}=s$ and $d\tilde{N}_{i,n+1}^{(k)}(s)=0$ otherwise.

It is shown in Zou (2015) that (8) can be rewritten as

[TABLE]

where \tilde{R}^{(k)}(s)=R^{(k)}(s)-\big{(}\inf_{t\geq 0}\{X^{(k)}(t)=0\}-s\big{)} denotes the difference between the $k$th patient’s residual life with transplantation at time $s$ and residual life if never transplanted. Note that the last term in (9) is invariant to allocation rules, and therefore the expected long-term average total life is equal to the long-term average of expected life gain from transplantation

[TABLE]

plus a constant that is invariant to allocation rules. For simplicity of notation, we will use (10) as the objective function in searching for optimal allocation rules.

3 Optimizing the Average Expected Life Gain

3.1 Formulation of the Optimization Problem

In studying the effect of allocation rule $\Gamma$ on the long-term average expected life gain, a pivotal quantity is the long-term average occupancy of health state $i$ for $i\in\mathcal{X}$:

[TABLE]

The long-term average of expected life gain can be expressed as a function of the long-term occupancy, and the occupancy can be shown to satisfy certain constraints - so that an upper bound for the long-term expected life gain may be found by maximizing over the long-term occupancy, subject to the constraints.

A proof is outlined in Zou (2015) for the proposition that, with any fair allocation rule, the limit in (11) exists so that $\pi_{i}(s)$ is well defined for all $i\in\mathcal{X}$ and $s\in[0,T]$. The existence of the limit follows from the fact that with a fair allocation rule, the waiting list process $\{W_{t}:t\geq 0\}$ is strong Markov with respect to filtration $\{\mathcal{F}_{t}\}$ and is positive recurrent. Therefore by the ergodic theory there exists a finite invariant measure $\sigma$ and $\pi_{i}(s)$ can be expresses in terms of $\sigma$.

For any fair allocation rule, there exists a corresponding long-term transplantation rate

[TABLE]

where $\mathbf{Q}=\{q_{ij}\}$ is the matrix of counterfactual rate of transition in health status. The transplantation rate is the difference between the rate of evolution that would occur absent organ allocation and the rate of the evolution of occupancy with allocation. For $s=0$, let $\pi_{i}(s-)=P(X(0)=i)$. For achievable long-term occupancies, $\pi_{i}(s)$ is smooth. The definition may be extended to the non-smooth case, which is relevant at the upper bound, by taking $\Psi_{i}(\{s\})=\pi_{i}(s-)-\pi_{i}(s)$ and allowing $\Psi_{i}(\{s\})>0$.

The effect of an allocation rule on the long-term average of expected life gain (10) can be expressed in terms of the long-term transplantation rate. Let $\tilde{\mu}_{i}(s)=E(\tilde{R}(s)\mid X(s)=i)$ and denote $\tilde{\boldsymbol{\upmu}}_{s}=(\tilde{\mu}_{1}(s),\dots,\tilde{\mu}_{n}(s))$. The result is stated formally in the following theorem.

Theorem 3.1**.**

[TABLE]

While each fair allocation rule has a corresponding transplantation rate $\boldsymbol{\Psi}$, not every $\boldsymbol{\Psi}$ can be traced back to an allocation rule. $\boldsymbol{\Psi}$ and $\boldsymbol{\pi}$ that correspond to fair allocation rules satisfy at least the following constraints. First, the total allocation rate has an upper bound associated with the rate of organ arrivals relative to the rate of patient arrivals:

[TABLE]

where $\mathbf{1}_{n}$ is the vector of length $n$ in which all elements equal to $1$. Essentially, the left side of (14) is the limiting proportion of patients who receive transplantation.

Second, let $\{\cdot\}_{i}$ denote the $i$th element of a vector, then for all $s\in[0,T]$ and $i\in\mathcal{X}$,

Proposition 3.2**.**

$\pi_{i}(s)\in[0,\pi_{i}(s-)]$, $d\pi_{i}(s)\in(-\infty,\{\boldsymbol{\pi}_{s-}\mathbf{Q}_{s}ds\}_{i}]$, 2. 2.

$\Psi_{i}(\{s\})\in[0,\pi_{i}(s-)]$, $d\Psi_{i}(s)\in[0,\infty)$, 3. 3.

$d\pi_{i}(s)\in[0,\{\boldsymbol{\pi}_{s-}\mathbf{Q}_{s}ds\}_{i}]$, $d\Psi_{i}(s)\in[0,\{\boldsymbol{\pi}_{s-}\mathbf{Q}_{s}ds\}_{i}]$, if $\pi_{i}(s-)=0$.

An upper bound for long-term average expected life is given by maximizing (13) with respect to $\boldsymbol{\Psi}$ and $\boldsymbol{\pi}$ subject to (14) and Proposition 3.2.

3.2 The Optimal Allocation Rate

This section focuses on the characterization of the optimal $\boldsymbol{\pi}$ and $\boldsymbol{\Psi}$ that maximizes the objective while satisfying the constraints. The primal-dual framework is used here to convert the constrained primal problem to an unconstrained dual problem, in which a penalty parameter $c$ associated with the constraint on the limiting proportion of transplanted patients is introduced. Then we apply the Hamilton-Jacobi-Bellman equations to the unconstrained problem and recursively solve for the optimal transplantation rate given any value of $c$. Dual-primal duality can be shown for this problem given the convexity of the objective and the constraints so that the primal problem is ultimately solved by finding the $c$ corresponding to the bound in the constraint.

Consider the unconstrained objective

[TABLE]

with a penalty parameter $c\geq 0$. Denote the $\boldsymbol{\Psi}$ that maximizes (15) while satisfying the boundedness constrains in Proposition 3.2 by $\boldsymbol{\Psi}^{c}$. The following result characterizes $\boldsymbol{\Psi}^{c}$ in a recursive manner.

Theorem 3.3**.**

For all $s\in[0,T]$ and $i\in\mathcal{X}$,

[TABLE]

if $\pi_{i}(s-)>0$, and

[TABLE]

if $\pi_{i}(s-)=0$. Here

[TABLE]

where $\eta_{i}^{c}(s+)$ is the $i$th element of the vector

[TABLE]

and $\gamma_{i}^{c}(s+)$ is the $i$th element of the vector

[TABLE]

where $d\boldsymbol{\Lambda}^{c}_{s}$ is a diagonal matrix with the $i$th diagonal element defined as $d\Lambda_{i}^{c}(s)=d\Psi^{c}_{i}(s)/\pi_{i}(s-)$ if $\pi_{i}(s-)\neq 0$ and $\Lambda^{c}_{i}(\{s\})=d\Psi^{c}_{i}(s)/\{\boldsymbol{\pi}_{s-}d\mathbf{Q}_{s}\}_{i}$ if $\pi_{i}(s-)=0$. The symbol $\prodi$ indicates the product integral operator as defined in Gill and Johansen (1990).

Remark*.*

Using the product integral $\prodi$ here enables a unified expression of the transition probability of the Markov process, whether $\boldsymbol{\Psi}$ is absolutely continuous or singular, or a mixture of both. For details on the application of $\prodi$ to the waiting list with allocations, see Zou (2015).

Here is a heuristic interpretation of Theorem 3.3. Imagine a scenario in which optimal allocation for a value $\rho$ of the organ arrival rate has been in place. The optimal allocation is characterized by $d\boldsymbol{\Lambda}^{c}$ where $c$ is associated with the allocation allowance $\rho$. Suppose we are now allowed an increase in the allocation rate to $\rho+d\rho$. It is shown in the Appendix that $\eta^{c}_{i}(s+)$ is in fact the conditional expected life gain for a patient in state $i$ at time $s$ past arrival, if not transplanted immediately but subject to future transplantation, and $\gamma^{c}_{i}(s+)$ is the corresponding probability of future allocation. If we choose to apply our $d\rho$ of allocation to subjects reaching state $i$ at time $s$, we will further free up $\gamma^{c}_{i}(s+)d\rho$ more allocation that could be applied at state $i$ and time $s$, further freeing up $(\gamma^{c}_{i}(s+))^{2}d\rho$, $\dots$ Ultimately, there would be a $d\rho/(1-\gamma^{c}_{i}(s+))$ increase in allocation at state $i$ and time $s$ and an expected life again of $\tilde{\mu}_{i}(s)-\eta^{c}_{i}(s+)$ per unit of increase in allocation. This suggests that optimal allocation occurs at those $i$ and $s$ such that

[TABLE]

is large.

4 Allocation Strategy

Theorem 3.3 suggests an upper bound of the long-term average of the expected life gain can be achieved by allocating organs to patients in state $i$ at time $s$ if and only if $\varphi_{i}^{c}(s)>c$, where $c$ is selected such that $\int_{[0,T]}d\boldsymbol{\Psi}^{c}_{s}\mathbf{1}_{n}=\rho/\tau$.

In reality, however, this upper bound cannot be reached, as when an organ is available for transplantation, there might not be any patient in the optimal state and waiting time in the current waiting list. Since the organ can only be preserved for a limited time before losing its functionality, the transplantation cannot be delayed for patients in the optimal state and waiting time to appear. One of the patients available in the waiting list, though in suboptimal states and waiting times, has to be selected for transplantation. Therefore, an allocation strategy is needed to prioritize patients in all possible states and waiting times.

Here we propose an allocation strategy motivated by the form of the optimal transplantation rate. The proposed strategy are developed based on three monotonicity results. First, the total transplantation rate is a monotone function of the penalty parameter $c$.

Proposition 4.1**.**

Suppose $c_{1}>c_{2}\geq 0$, then

[TABLE]

Second, the long-term average expected life gain from transplantation is a monotone function of $c$.

Proposition 4.2**.**

Suppose $c_{1}>c_{2}\geq 0$, then

[TABLE]

Third, the optimal set of state and waiting time pairs for transplantation is monotone with respect to $c$. Formally, for each value of $c\geq 0$, define the corresponding allocable set

[TABLE]

which consists of state and waiting time pairs with non-zero transplantation rate in the optimal solution given by Theorem 3.3. The following result states that $\mathcal{A}_{c}$ expands monotonically as $c$ decreases.

Theorem 4.3**.**

Suppose $c_{1}>c_{2}\geq 0$, then

[TABLE]

Based on the monotonicity results, a full order priority ranking of all state and waiting time pairs can be obtained with the following algorithm outlined in the pseudo code Algorithm 1 below: start with a large enough value of $c$ such that $\mathcal{A}_{c}=\emptyset$, then the value of $c$ is gradually decreased, representing increasing organ availability, and the order in which each state and waiting time pair enter the allocable set $\mathcal{A}_{c}$ as $c$ decreases is recorded. A binary search is used to find each value of $c$ when a new pair of $(i,s)$ enters $\mathcal{A}_{c}$. Pairs that enter the allocable set earlier are given higher priorities for transplantation, as they are selected for allocation to prolong the long-term average expected total life when the availability of organs is more stringent and remain allocable when more organs are available. The sensitivity parameter $\alpha$ in Algorithm 1 is a small positive number needed in the algorithm to avoid infinite loops. Smaller $\alpha$ leads to higher precision in ranking all waiting time and state combinations. $N_{\text{steps}}$ is the maximal number of steps allowed in calculating the ranks of $(i,s)$ combinations.

Once the rank of state and waiting time pairs is determined, whenever an organ is retrieved, patients available in the current waiting list can be instantly ranked based on their state and waiting time. The patient with the highest rank will be selected for transplantation. The resulting long-term average expected total life using this strategy is expected to be close to the upper bound, as patients are ranked according to their potential contributions to the limiting average total life considering the limited availability of organs. This allocation strategy satisfies the fairness requirement, as the allocation decisions are solely determined by states and waiting times of patients in the current list.

5 Matching of Multiple Organ Types

Now we generalize the allocation problem with single organ type to multiple organ types. Suppose there are $N_{o}>1$ types of organs with corresponding arrival rates $\{\rho^{(w)}:w=1,\dots,N_{o}\}$. The interesting case is when $\sum_{w=1}^{N_{o}}\rho^{(w)}<\tau$ which reflects the reality that there are not enough organs to meet the need of the patients.

For any $w\in\{1,\dots,N_{o}\}$, let $\tilde{\upmu}_{i}^{(w)}(s)$ denote the expected gain in residual life for a patient in state $i$ and waiting time $s$ with immediate transplantation with a type $w$ organ comparing to never transplanted. Note $\tilde{\upmu}_{i}^{(w)}(s)$ is a function of both the state and waiting time, providing the flexibility to model the benefit of transplantation with each organ type as depending on the compatibility of the organ with the patient and possibly varying with waiting time. In addition to factors that are usually taken into account in matching donors and recipients such as blood type, patients’ health states can include factors that are found to be significantly related to the benefits of different organ types. For example, studies (Oto et al. (2004), Bonser et al. (2012), Shigemura et al. (2013) and Kopp et al. (2015)) have found that while the post-transplant survival is worse with organs from donors with a smoking history, certain recipient characteristics could make a difference in the benefit of transplantation with a smoking-donor organ.

5.1 The Optimal Transplantation Rate for Multiple Organ Types

It can be shown, as in the single organ type case, that with any allocation rule that satisfies fairness requirements for each organ type there exists a limiting transplantation rate $\boldsymbol{\Psi}^{(w)}_{s}$ for each $w=1,\dots,N_{o}$, with which the limiting averaged life gain can be represented by

[TABLE]

Moreover, for each $w$, the constraint

[TABLE]

must hold as the limiting proportion of patients transplanted with organ type $w$ is bounded by the ratio of type $w$ organ arrival rate to the patient arrival rate. Each $\boldsymbol{\Psi}_{s}^{(w)}$ must also satisfy the first condition in Proposition 3.2 and the following constraints:

[TABLE]

and

[TABLE]

as a patient can only be transplanted with one organ type at a time.

As in the single organ type case, the optimal allocation rates can be solved by first considering the penalized objective

[TABLE]

where each $c^{(w)}\geq 0$ is a penalty parameter associated with the constraint (25), indicating the availability of type $w$ organs. Let $\mathbf{c}=(c^{(1)},\dots,c^{(N_{o})})$ denote the vector of penalty parameters for all organ types. The following result characterizes the optimal allocation rates for (28) given any fixed value of $\mathbf{c}$.

Theorem 5.1**.**

For all $s\in[0,T]$ and $i\in\mathcal{X}$, the transplantation rates $\{\Psi_{i}^{(w),\,c}:w=1,\dots,N_{o}\}$ that maximize (28) are given by

[TABLE]

if $\pi_{i}(s-)>0$, and

[TABLE]

if $\pi_{i}(s-)=0$, where $\eta_{i}^{(w),\,c}(s+)$ is the $i$th element of the vector

[TABLE]

and $\gamma_{i}^{(w),\,c}(s+)$ is the $i$th element of the vector

[TABLE]

where each $d\boldsymbol{\Lambda}^{(w),\,c}_{s}$ is a diagonal matrix with the $i$th diagonal element being $d\Lambda^{(w),\,c}_{i}(s)=d\Psi^{(w),\,c}_{i}(s)/\pi_{i}(s-)$ if $\pi_{i}(s-)\neq 0$ and $\Lambda^{(w),\,c}_{i}(\{s\})=d\Psi^{(w),\,c}_{i}(s)/\{\boldsymbol{\pi}_{s-}d\mathbf{Q}_{s}\}_{i}$ if $\pi_{i}(s-)=0$.

Theorem 5.1 indicates that the optimal allocation rates for multiple organ types in the heavy traffic setting can be achieved by allocating all of patients’ occupancy in state $i$ at waiting time $s$ to the transplanted state with type $w$ organs if both of the following conditions are satisfied: first, the difference between the post-transplant residual life gain and the organ-type-specific penalty, if immediately transplanted, is the largest when transplanted with a type $w$ organ comparing to other organ types, which makes type $w$ the “designated” organ type for $(i,s)$; and second, the difference between the post-transplant residual life gain and the organ-type-specific penalty parameter, if transplanted with a type $w$ organ immediately, is larger than the expected difference between the residual life gain and the penalty parameter if the patient is not transplanted immediately but open to future transplantation opportunities with any type of organ.

5.2 Allocation Strategy for Multiple Organ Types

In practice, as in the single organ type case, the optimal transplantation rates usually cannot be reached, as when a type $w$ organ becomes available, there may not be any patient in the optimal state and waiting time for this specific organ type in the waiting list. Consequently, an allocation strategy is needed to rank patients’ priorities of transplantation with each organ type.

Theorem 5.1 motivates a natural allocation strategy. Given the form of the optimal transplantation rate, for each $w\in\{1,\dots,N_{o}\}$, the allocable set of state and waiting time combinations for type $w$ organs is denoted by

[TABLE]

in which $\mathbf{c}$ indicates the availability of each type of organs. If each $\mathcal{A}^{(w)}_{\mathbf{c}}$ expands monotonically when $\mathbf{c}$ decreases to reflect more organs becoming available, then a priority ranking for transplantation with each organ type can be obtained by comparing the order of each $(i,s)$ entering the $\{\mathcal{A}^{(w)}_{\mathbf{c}}\}$ as $\mathbf{c}$ decreases. Pairs that enter $\mathcal{A}^{(w)}_{\mathbf{c}}$ earlier are given higher priorities for transplantation with type $w$ organ.

The monotonicity of $\{\mathcal{A}^{(w)}_{\mathbf{c}}\}$, however, generally does not hold when $\mathbf{c}$ decreases freely. In fact, the monotonicity in the multiple organ type case only holds in a weaker sense, that is, when the decreasing of $\mathbf{c}$ follows special “paths” with constraints on changes in each $c^{(w)}$. The following results describe two types of weak monotonicity properties. First, when all of the $\{c^{(w)}\}$ decrease by the same amount, the allocable set for each organ type expands monotonically. Second, when fixing all the $c^{(v)}$ where $v\neq w$ and decreasing only $c^{(w)}$, the allocable set for type $w$ organs expands monotonically.

Proposition 5.2**.**

If $c^{(w)}_{1}-c^{(w)}_{2}$ is a positive constant independent of $w$, then $\mathcal{A}^{(w)}_{\mathbf{c}_{1}}\subseteq\mathcal{A}^{(w)}_{\mathbf{c}_{2}}$ for all $w\in\{1,\dots,N_{o}\}$.

Proposition 5.3**.**

If $c^{(w)}_{1}>c^{(w)}_{2}$ and $c^{(v)}_{1}=c^{(v)}_{2}$ for all $v\neq w$, then $\mathcal{A}^{(w)}_{\mathbf{c}_{1}}\subseteq\mathcal{A}^{(w)}_{\mathbf{c}_{2}}$.

Following the above paths of decreasing $\mathbf{c}$, the priority rankings for all organ types can be obtained using a three-stage algorithm:

The value of $\mathbf{c}$ such that the active constraints (25) are satisfied, at least approximately so, are calculated via a multiple variable gradient descent algorithm. This step guarantees that $\mathbf{c}$ reflects the practical relative availability of each organ type. 2. 2.

Given the calculated value of $\mathbf{c}$, increase $\mathbf{c}$ by a constant so that the new $\mathcal{A}^{(w)}_{\mathbf{c}}$ is empty for each $w$. Then decrease $\mathbf{c}$ by the same amount as in Proposition 5.2 and rank the priorities for type $w$ organ transplantation according to the order of each $(i,s)$ pair entering $\mathcal{A}^{(w)}_{\mathbf{c}}$.

When decreasing each $c^{(w)}$ by the same amount, the designated organ type $w$ for each $(i,s)$ that maximizes $\tilde{\upmu}^{(w)}_{i}(s)-c^{(w)}$ remains the same. Consequently, when $\mathbf{c}$ is small enough for each $(i,s)$ to be allocable, each $(i,s)$ has a designated organ type and the space of $\mathcal{X}\times[0,T]$ is partitioned into disjoint subsets indexed by the designated organ types. 3. 3.

For each $w$, decrease $c^{(w)}$ further while keeping all other $c^{(v)},v\neq w$ constant as in Proposition 5.3, until all $(i,s)\in\mathcal{X}\times[0,T]$ are ranked for type $w$ organs. This step guarantees that for each organ type there is a complete ranking of combinations in $\mathcal{X}\times[0,T]$. This step is necessary since in practice when an organ of type $w$ arrives, there might not be any patient in the subset of $\mathcal{X}\times[0,T]$ with type $w$ as the designated organ type. An allocation decision still needs to be made, since there cannot be a delay in transplantation due to the time limit for the organ to keep its functionality, and patients who are currently in the waiting list might still benefit from the transplantation. Complete rankings of $\mathcal{X}\times[0,T]$ for each organ type ensures an allocation decision can be made to select the patient who would benefit the most when an organ becomes available.

6 Comparison of Allocation Strategies

This section focuses on the application of the proposed allocation strategy in the single-organ-type scenario in the context of realistic models for health state transitions of patients in the lung transplantation waiting list. Waiting lists are simulated based on parameters estimated from the United Network for Organ Sharing (UNOS) lung transplantation data, and the proposed strategy is compared to the Lung Allocation Score (LAS) system by its performance in extending life when applied to the simulated waiting lists.

6.1 The United Network for Organ Sharing (UNOS) Data

The data used are the waiting list, transplant and follow-up UNOS Standard Transplant Analysis and Research files, which contain information of heart, lung, and simultaneous heart-lung registrations and transplants that were listed or performed in the United States and reported to the Organ Procurement and Transplantation Network (OPTN) from October 1, 1987 to December 31, 2012. In the data, the first transplantation with the LAS occurred on May 5, 2005. Before June 11, 2013, only patients at least 12 years of age received priority for deceased donor lung offers based on the LAS, and thus we only include patients at least 12 years of age in the data analysis. There were 16,049 such patients registered in the above time period, with 129,881 records of medical measurements updated sporadically at different times for different patients during their tenure on the waiting list. Among the 16,049 candidates, 64.6% received transplantation, 18.4% died while waiting in the list, and the remaining were still waiting at the end of the study period or censored by loss of follow-up. Among patients who received transplantation, 37.8% died and 62.2% were still alive at the end of the study period or were censored due to loss of follow-up.

6.2 Parameter Estimation from the UNOS Data

We fit two separate proportional hazards regression models with time varying covariates to estimate the time-on-waiting-list specific hazards for in-waiting-list and post-transplant deaths. This serves two purposes: 1. to calculate the LAS given a patient’s current covariates values and 2. to characterize patients’ health states.

Covariates in the proportional hazard regressions are those used in the calculations in the revised LAS system by UNOS (for details of covariates used in the LAS calculation, see UNOS (2017)), whenever they are available in the data, so that no advantage is gained from an improved variable selection when comparing the proposed strategy to the LAS. Note in estimating the hazard for in-waiting-list death, transplantation leads to censoring in the data. The censoring by transplantation can be treated as censoring at random, nevertheless, as the selection of organ recipients was based on patients’ LAS calculated with the same covariates used in estimating the hazard for in-waiting-list death and thus can be treated as conditionally independent of future survivals given the observed covariates. See Tables LABEL:t1 and LABEL:t2 in the Appendix for lists of covariates and coefficient estimates in the two proportional hazards regressions. More details including estimated baseline survival curves can be found in the Supplementary Material.

6.3 Defining Patients’ Health States

Patients’ health states are characterized in terms of the linear combinations of the covariates in the proportional hazards models, with the estimated regression coefficients being the weights associated with the covariates. Each patients’ health state at waiting time $s$ is defined as the pair $(S_{wl}(s),S_{\mu}(s))$, where $S_{wl}(s)$ and $S_{\mu}(s)$ are, respectively, values of the linear combinations of covariates at $s$ in the proportional hazards models for in-waiting-list and post-transplant survivals.

Characterizations of transitions in health states are based on observations from exploratory analysis of the data. Covariates relevant to waiting list and post-transplant survivals are categorized into deterministic and stochastic variables. Deterministic variables are those remain mostly constant or are deterministic functions of waiting time. Covariates such as age, diagnosis group and detailed diagnosis fall into this category. Stochastic variables are those change randomly while patients are waiting for transplants. Examples of stochastic variables include BMI, functional assistance status, ventilation status, creatinine, oxygen and six minute walk distance.

In the LAS calculation (UNOS (2017)), the set of the stochastic covariates for estimating post-transplant survival is a subset of the stochastic covariates used in estimating the in-waiting-list survival. Therefore, the set of stochastic covariates can be partitioned further into two sub-categories: covariates used in both proportional hazards models for in-waiting-list survival and post-transplant survival, and covariates used only for in-waiting-list survival. As a result, the linear combinations $S_{wl}(s)$ and $S_{\mu}(s)$ can be written as

[TABLE]

where $\mathbf{X}_{1}$ represents stochastic covariates included in both models, $\mathbf{X}_{2}$ represents stochastic covariates for $S_{wl}$ only, $\mathbf{X}_{3}$ represents deterministic covariates for $S_{wl}$, and $\mathbf{X}_{4}$ represents deterministic covariates for $S_{\mu}$. See Tables LABEL:t1 and LABEL:t2 in the Appendix for detailed categorizations of the covariates into each of the five $\mathbf{X}$ categories.

$\boldsymbol{\upbeta}_{1}\mathbf{X}_{1}(s)$ and $\boldsymbol{\upbeta}_{2}\mathbf{X}_{2}(s)$ are characterized as following conditionally independent compound jumping processes whose jumping intensities and magnitudes depend on the current waiting time and health state $(S_{wl},S_{\mu})$. The jumping hazards are estimated with proportional hazards regressions for recurrent events, while the jumping magnitudes (in logarithmic scale) are estimated with linear regressions. It was observed from the data that estimates of corresponding coefficients in $\boldsymbol{\upbeta}_{1}$ and $\tilde{\boldsymbol{\upbeta}}_{1}$ are approximately proportional (Figure 1). Therefore the transition in $\tilde{\boldsymbol{\upbeta}}_{1}\mathbf{X}_{1}(s)$ can be approximated by the transition in $\boldsymbol{\upbeta}_{1}\mathbf{X}_{1}(s)$ scaled by the estimate of their ratio, and vice versa.

The transition probabilities in $(S_{wl}(s),S_{\mu}(s))$ can be derived from the transition probabilities of $\boldsymbol{\upbeta}_{1}\mathbf{X}_{1}(s)$ and $\boldsymbol{\upbeta}_{2}\mathbf{X}_{2}(s)$ and the deterministic trajectories of $\boldsymbol{\upbeta}_{3}\mathbf{X}_{3}(s)$ and $\boldsymbol{\upbeta}_{4}\mathbf{X}_{4}(s)$. Specifically, for any $U,V\subset\mathbb{R}$,

[TABLE]

where $V_{x}$ denotes the interval $V$ shifted by $-x$ and $l(\Delta\tilde{\boldsymbol{\upbeta}}_{1}\mathbf{X}_{1}(s))=\Delta\boldsymbol{\upbeta}_{1}\mathbf{X}_{1}(s)$ is the linear function approximating the linear relationship between the changes in $\boldsymbol{\upbeta}_{1}\mathbf{X}(s)$ and $\tilde{\boldsymbol{\upbeta}}_{1}\mathbf{X}(s)$, and each conditional probability in the right side of (34) is calculated using the estimated jumping intensities and magnitudes described earlier. The distribution of initial states $(S_{wl}(0),S_{\mu}(0))$ is estimated with the empirical distribution of patients’ states upon arrival to the waiting list.

6.4 Discretization of Patients’ Health States

For the purpose of efficiency in computation, we discretize waiting time to periods of $30$ days and $(S_{wl},S_{\mu})$ to a $4\times 4$ finite state space. Each of $S_{wl}$ and $S_{\mu}$, which were originally in continuous scales, are partitioned into $4$ intervals, where the first and last cut-off points of the intervals are the 1/4 and 3/4 empirical quantiles of the original variables defined in (33), and the rest of the cut-off points are such that the intervals have equal lengths. As a result there are $4$ states: $\{1,2,3,4\}$ for each of $S_{wl}$ and $S_{\mu}$, and the overall discretized patient state space is defined as the cross product of the discretized state spaces of $S_{wl}$ and $S_{\mu}$.

The corresponding discrete transition probabilities among states are calculated from the estimated continuous transition probabilities. Specifically, at each waiting time $s$, the transition probability from state $i=(i_{1},i_{2})$ to $j=(j_{1},j_{2})$ as waiting time increases by $\Delta s$ is approximated with (34), in which $x$ and $y$ are the middle points (medians in case of the first and last intervals) of the $i_{1}$th interval of $S_{wl}$ and $i_{2}$th interval of $S_{\mu}$, and $U$ and $V$ are the $j_{1}$th interval of $S_{wl}$ and $j_{2}$th interval of $S_{\mu}$, respectively. We first estimate the single-day transition probabilities with $\Delta s=1$, and then the cumulative transition probabilities for the discretized waiting time intervals with $\Delta s=30$ are calculated as the product of the single-day transition probabilities using the Markov property.

We also experimented with different fineness levels in the discretization of waiting time and health states as well as different models for transitions in health states and simulation results were robust to different settings. See the Supplementary Material for another example from Zou (2015) in which waiting time was discretized into $90$-day periods and $(S_{wl},S_{\mu})$ was discretized into a $3\times 3$ state space and the transition probabilities of the discretized $(S_{wl},S_{\mu})$ was estimated with a multinomial regression model.

6.5 Calculation of Priority Rankings of Different Allocation Strategies

Given the defined health states, allocation strategies can be formulated as priority rankings of all combinations of waiting time and health state. Here we compare our proposed allocation strategy with the LAS and the refined LAS-type method. First, the LAS currently used by the UNOS is $\textrm{LAS}=100\cdot(\textrm{PTAUC}-2\cdot\textrm{WLAUC}+730)/1095$, where WLAUC is the estimated in-waiting-list life expectancy during an additional year and PTAUC is the estimated post-transplant life expectancy during the first year, given the patient’s current waiting time and state $(S_{wl},S_{\mu})$. Patients with higher LAS have more priorities of transplantation. Specifically, the two life expectancy measures are calculate as follows:

[TABLE]

where $\bar{F}_{wl}(t)$ is the in-waiting-list survival function at time $t$ (treating current waiting time as time [math]), and $\bar{F}_{tx}(t)$ is the post-transplant survival function at time $t$ (treating time at transplant as time [math]).

It is observed from the data that patients’ in-waiting-list and post-transplant life are usually much longer than one year. Studies showed the emphasis on one-year survival by the LAS might have led to worse long-term survival. See for example, Maxwell et al. (2014), for detailed statistics. Moreover, it is implicitly assumed in the LAS calculation that the in-waiting-list survival functional $\bar{F}_{wl}(t)$ is invariant to the patient’s current waiting time, which may not be case. Here we also calculate a refined LAS without the one-year constraint on the life expectancy measures and without assuming the invariance to current waiting time. Due to heavy censoring of organ recipients with post-transplant residual life longer than five years, in the refined LAS we calculated the median of the covariate-specific post-transplant survival, as it is less sensitive to missing values comparing to the mean, as the post-transplant life measure.

The priority ranking of all waiting time and health state combinations in our proposed allocation strategy is calculated with the following procedure. Given parameters estimated from the UNOS data, for a fixed value of the penalty parameter $c$, whether $(s,i)$, in which $s$ is the waiting time and $i$ is the value of the health state $(S_{wl},S_{\mu})$, is allocable is decided by whether the optimal transplantation rate for this combination is non-zero. The optimal transplantation rate can be calculated with formula (17), which in turn is decided by $\varphi_{i}^{c}(s)$ in (18). The calculation of $\varphi_{i}^{c}(s)$ requires the estimation of $\tilde{\mu}_{i}(s)$, $\eta_{i}^{c}(s+)$, and $\gamma_{i}^{c}(s+)$. $\tilde{\mu}_{i}(s)$ is calculated as the difference between the estimated post-transplant residual life and the estimated in-waiting-list life if never transplanted, where both residual life are estimated with proportional hazards models given the current waiting time and health state. The estimated $\eta_{i}^{c}(s+)$ and $\gamma_{i}^{c}(s+)$ are calculated using formulae (19) and (20), in which the counterfactual transition matrix $\mathbf{Q}$ is estimated from the data with the method described in Sections 6.3 and 6.4 and $d\mathbf{\Lambda}^{c}$ is calculated backwards in an iterative manner from the upper bound $T$ of waiting time. Finally, based on the monotonicity of the allocable sets $\mathcal{A}_{c}$ in Theorem 4.3, a full order ranking of transplant priorities for combinations of waiting times and states is obtained by comparing the order of their entries to $\mathcal{A}_{c}$ when $c$ decreases as described in Algorithm 1.

See Figures 2, 3 and 4 below for a comparison of priority rankings with the proposed allocation strategy and the two LAS-type methods. In each of the graphs, numbers $1-100$ floating on the gray background indicate the discretized waiting time periods of $30$ days. Each mosaic pattern under the waiting time title contains the $4\times 4$ state space of $(S_{wl},S_{\mu})$, where the $x$-axis indicates the states of $S_{wl}$ and the $y$-axis indicates the states of $S_{\mu}$. Larger index of the states corresponds to shorter life. Therefore, each small rectangle in the graph shows the priority of a combination of waiting time and health states, where blue indicates higher priority for transplantation, while red indicates lower priority. For example, the rectangle on the very top-left of Figure 2 is in deep red, demonstrating that with the proposed allocation strategy, patients that are in the first $30$-day waiting time period and with $(S_{wl},S_{\mu})=(1,4)$ have very low priority of transplantation. This is reasonable, since $S_{wl}=1$ indicates a prediction of long in-waiting-list residual life while $S_{\mu}=4$ predicts a short post-transplant residual life, and thus low benefit of transplantation, especially given the fact that these patients are new to the waiting list and are less urgent in receiving transplants.

6.6 Comparison of Policy Performances via Simulation Results

Finally, we simulate waiting lists with parameters estimated from the UNOS data. Specifically, patient and organ arrivals are simulated according to independent homogeneous Poisson processes. The organ arrival rate $\rho$ and the patient arrival rate $\tau$ are estimated with average numbers of patient arrivals and organ arrivals per $30$-day period after year 2006: $\hat{\rho}=104$, $\hat{\tau}=173$. At the end of each waiting time period, counterfactual transitions in patients’ health states are simulated according to the estimated transition probabilities. Using priority rankings of the proposed allocation strategy and the two LAS methods, organs, whenever available, are allocated to patients with the highest ranked combinations of waiting times and states in the current waiting list. Patients who are not selected for transplants remain on the waiting list with possibilities of future counterfactual transitions and transplantations. At the end of the simulation, the total life, life in the waiting list and post-transplant life, averaged over all patients ever entered the waiting list, are counted in days for each of the allocation strategies.

Table 1 reports the means and standard deviations of averaged life outcomes of $200$ independent simulations. In each of the $200$ simulation runs, a waiting list containing $1605$ patients was generated and $973$ organs were allocated to the patients with each allocation strategy. In addition to the above mentioned strategies, we consider a random allocation strategy in which patients are chosen for transplantation randomly. In each simulation run, the generated waiting list is copied four times to apply each of the allocation strategies in question. The averaged total life under each allocation strategy in Table 1 is also compared to the model upper bound of $1881$ days.

By applying the proposed allocation strategy, the averaged total life increased by $7.7\%$ comparing to the result with the LAS. The gain in total life came from a much improved post-transplant life: it increased by $16.8\%$ comparing to the post-transplant life under the LAS. Meanwhile, the average time in waiting list for patients (including those who never received transplantation and the organ recipients) was shortened by $11.7\%$. Applying the refined LAS without the one-year constraint in the LAS calculation and without the assumption that the waiting-list survival is invariant to the current waiting time slightly improves the averaged total life. As expected, the purely random allocation had a much worse performance than the proposed strategy and the LAS methods.

Figure 5 shows the box plots of the averaged life outcomes with each allocation strategy. From left to right the averaged total life with the model upper bound of $1881$ (red dashed line), the averaged life in waiting list, and the averaged life post transplant are displayed. An interesting observation is that while the randomized strategy resulted in a larger variance, especially in post-transplant life, which is expected given that organ recipients were selected randomly, the mean post-transplant life of patients under the random strategy was not substantially different from those under the LAS methods. This observation might indicate the current LAS system still has much space to improve, especially in extending patients’ post-transplant life.

Regarding patients’ average life in the waiting list with each allocation strategy, Table 2 and Figure 6 show detailed average waiting time for organ recipients and patients who never received transplantation.

For patients who died while waiting for transplantation, the original LAS with one-year constraint led to the longest in-waiting-list life, possibly due to its emphasis on in-waiting-list survival. The LAS gives higher priority of transplantation to patients with large probability of dying within one year. Without the one-year constraint, the refined LAS allocated less organs to those who are likely to die in one year without transplant and hence the shorter average in-waiting-list life. Overall, the LAS methods emphasize more on the in-waiting-list survival over the post-transplant survival.

On the other hand, the proposed strategy aims to optimize patients’ total life and thus patients with larger expected life extension with transplantation were given higher priority. Since some of the patients who were likely to die in one year without transplantation were also expected to have short post-transplant life, they were not selected as organ recipients, which resulted in a shorter average waiting list life.

For patients who received transplantation, the proposed method shortened the average waiting time for an organ comparing to the LAS methods, which might be considered to be an advantage by the clinicians and patients. The waiting time was shorter with the original LAS comparing to the refined LAS, which was possibly due to allocations of organs at early waiting time to patients who were otherwise going to die soon.

7 Numerical Results on Multiple Organ Type Matching

Obtaining substantive results in the multiple organ type scenario would require a comprehensive study of organ properties from real transplantation data. Here we demonstrate the proposed allocation strategy for multiple organ types with a synthetic data example.

Two organ types are generated in the example. The post-transplant residual life with type-$1$-organ transplantation is larger for all $(i,s)$ combinations comparing to with type-$2$-organ transplantation. However, for some $(i,s)$ combinations, the difference in the benefit between the two organ types is less significant, modeling the situations where patients in certain states at certain stage of the disease may have a better compatibility with type $2$ organs. This example is generated to resemble the phenomenon that patients generally have a better survival if transplanted with non-smoking-donor organs, but recipients with certain characteristics may suffer less from the reduction in benefit when transplanted with smoking-donor organs.

Following the three-stage algorithm in Section 5.2 motivated by the monotonicity properties Propositions 5.2 and 5.3, the partition of the space of $(i,s)$ combinations into disjoint subsets by their designated organ types is obtained, and complete rankings of priorities for each organ type are calculated. All results are calculated with $\sum_{w=1}^{2}\rho^{(w)}/\tau=104/173$ and $\rho^{(1)}/\rho^{(2)}=1$.

Results are visualized in Figures 7 - 10. In each of the graphs, numbers $1-100$ floating on the gray background indicate the discretized waiting time periods of $30$ days. Each mosaic pattern under the waiting time title contains the $3\times 3$ state space of $(S_{wl},S_{\mu})$, where the $x$-axis indicates the states of $S_{wl}$ and the $y$-axis indicates the states of $S_{\mu}$.

Figure 7 shows the partition of $(i,s)$ combinations by designated organ type. The majority of $(i,s)$ combinations have color blue, indicating type 1 (non-smoking-donor) organ is the designated organ type. However, for some $(i,s)$ combinations the designated organ type is type $2$ (smoking-donor) organ. This partition comports with the differences between post-transplant residual life with type $1$ and type $2$ organs visualized in Figure 8. In particular, positions of the white rectangles in Figure 7 match those of the dark gray rectangles in Figure 8, which are $(i,s)$ combinations that suffer less from the reduction in post transplant residual life with type $2$ organs.

Figures 9 and 10 visualize the priority rankings of $(i,s)$ combinations for transplantation with type $1$ and type $2$ organs respectively. In Figure 9 it is expected that the priority is lowest whenever the color is white in the partition in Figure 7, as those $(i,s)$ combinations have designated organ type 2 and will only be considered for type-$1$-organ transplantation if no patients in $(i,s)$ colored blue, indicating type 1 as designated organ type, are available in the current waiting list when a type 1 organ arrives. Among those $(i,s)$ combinations with designated organ type 1, the priority ranking of transplantation generally resembles the ranking in the single-organ case: $(i,s)$ with longer estimated post transplant residual life and shorter in-waiting-list residual life, such as the bottom-right rectangle in each mosaic, generally has higher priority of transplantation. Similar phenomena can be observed in the priority ranking of type 2 organ transplantation in Figure 10.

However, the priority rankings in the multiple organ type scenario do not always follow those in the single organ type case. “Flips” of the priorities are observed: some $(i,s)$ combinations with high priorities in the single organ type case now have lower priorities comparing to other combinations. Examples of such flips include the bottom row of mosaic No. 83 and mosaics No. 87 and 88 in Figure 9. This phenomenon of flips echos our discussion on the monotonicity in Section 5.2: the priority ranking for one organ type is affected by the relative availability of all organ types represented by $\mathbf{c}$ and the strong monotonicity of priorities usually does not hold when $\mathbf{c}$ decreases freely.

Overall, the above results demonstrate that multiple organ type allocation and matching is not a trivial extension of the single organ type problem. The ranking of $(i,s)$ combinations for each organ type depends on the partition of $(i,s)$ combinations by their designated organ types. Moreover, the partition and priority rankings are determined by the complex interplay between transitions in health states, the residual life of each state at each waiting time, and the relative availability of different organ types.

8 Concluding Remarks

We have presented an approach to modeling the lung transplantation waiting list and comparing allocation rules. Here patient and organ arrivals are modeled as independent homogeneous Poisson processes, and counterfactual patient health status trajectories absent transplantation are modeled as independent and identically distributed inhomogeneous Markov processes. Patients’ expected post-transplantation residual life depends on both the health state at the time of transplantation and the waiting time at transplantation. The model setting here is capable of capturing the randomness in patient and organ arrivals and the complex dynamics of patients’ health characteristics and their effects on life outcomes. In practice, the researchers and policy makers can always expand the state space of patient health state so that trajectories of the transitions of health states approximately follow the Markov property. For example, previous medical records of patients can be included if they are believed to contain important information for predicting future transitions.

Allocation rules are modeled as index sequences of transplanted patients. Only fair allocation rules are considered in the comparison. Under fair allocation rules, the choice of patients for transplantations are decided by patients’ health states and waiting times at the time of organ arrival and a random variable that is conditionally independent of patients’ past and future states and survivals. Allocation probabilities are also required to be invariant to patient index and calendar time. This definition of fairness has two implications: first, unrealistic rules that can predict patients’ future states and survivals are not considered; second, given patients’ health states and waiting times, the allocation decisions are independent of other factors. Therefore the definition of fairness here implicitly addresses the equity issue in organ allocations.

It is shown each fair allocation rule has a corresponding limiting transplantation rate. Under a fair allocation rule, the average rate of transplantation to patients in any state and waiting time converges to the corresponding transplantation rate as calendar time increases. The limiting transplantation rate satisfies constraints that reflect the scarcity of organs. The main constraint is that the average proportion of transplanted patients in the limit is bounded by the ratio of organ arrivals to that of patient arrivals. The limiting average total life (or equivalently, the limiting average life gain from transplantation), represented in terms of the transplantation rate, is used as the standard in comparing fair allocation rules.

The optimal transplantation rate subject to the constraints is characterized recursively with the Hamilton-Jacobi-Bellman equations. Then a fair allocation strategy is developed based on the form and monotonicity properties of the optimal transplantation rate. The allocation strategy is to use the penalty parameter $c$ associated with the constraint on the average proportion of transplanted patients as an index to prioritize patients’ states and waiting times.

The index $c$ is partially related to the Gittins Index (Gittins (1989)). Bertsimas and Niño-Mora (2000) considered a setting that is close to the organ transplantation problem, in which a fixed number of resources are to be allocated to a fixed number of subjects that are in constant transitions of states in discrete time. Allocating resources using Gittins Index leads to optimal or asymptotically optimal objectives in problems where subjects remain static if not selected for allocations (Weber (1992), Whittle (1980), Whittle (1988)) and in settings where the rates of patient and organ arrivals tend to infinity. In restless bandits problems in which all subjects are constantly in transitions, including the problem studied here, though there are sufficient conditions for Gittins Index to be optimal (Bertsimas and Niño-Mora (1996), Niño-Mora (2001)), the optimality is not guaranteed in general.

Existing restless bandits models do not cover all complexities of the lung transplantation allocation problem. These complexities include the random arrival of patients and organs to the waiting list, random exits of patients due to death and transplantation, constraints on the allocation when an organ becomes available due to limited availability of patients in allocable states, and fairness considerations.

In our paper we model the lung transplant problem accounting for these complexities and study asymptotic properties of the waiting list. Then we provide a heavy traffic solution to optimize patients’ long-term average life. The heavy traffic optimal solution given in Section 3 can be seen as a much generalized version of the Gittins Index: instead of having a discrete number of jobs in the restless bandits problem, we have a continuum of “jobs” (expected occupancy of states at waiting times) and the resources are allocated to the continuum of jobs under constraints that are also continuous. As discussed in Section 4, the heavy traffic optimal solution cannot be directly applied to the allocation of organs in reality, and we propose an allocation strategy motivated by the heavy traffic solution that can lead to a long-term average life that is close to the heavy traffic upper bound.

In short, part of our development is tangentially related to the restless bandits problem and the Gittins Index. The scope of the lung transplantation problem itself, however, is beyond the idealized restless bandits model due to its complexities.

Simulation studies show it may be possible to improve the Lung Allocation Score (LAS) currently used by UNOS and increase the average total life by as much as 7.7%. Results provided here are provisional and a deeper understanding of the lung allocation procedure and the optimal allocation strategy requires further effort. As discussed previously, there may be a gap between the objective using the proposed strategy and the practical upper bound of the objective. The gap may be a result of the lack of patients in optimal states and waiting times when organs are available, which may stem from the fact the constraints in the optimization problem do not cover all practical confinements on the allocations.

In future studies we will extend the donor-recipient matching results developed here and conduct comprehensive substantive studies using not only patients but also donor records from UNOS. Issues related to cross-region transplantation, or other practical aspects in lung transplantation may also require further constraints to be imposed in formulating the optimal transplantation rate and developing allocation strategies.

Acknowledgement

The authors would like to thank Eric Peterson for his assistance with the UNOS data and Marco Scarcini for discussions on this problem. We would also like to thank the AE and the reviewers for their invaluable comments and suggestions in the process of improving this paper.

Appendix

8.1 Coefficient Estimates for In-Waiting-List and Post-Transplant Survivals

8.2 Proofs of Main Technical Results

Here we outline the proofs of selected main results. For detailed proofs of all technical results, see Zou (2015).

Proof of Theorem 3.3

Proof.

For any $i\in\mathcal{X}$, $\Psi_{i}$ can be decomposed into (c.f. for example Halmos (1974)),

[TABLE]

where $\Psi_{i}^{ac}$ is absolutely continuous with respect to the Lebesgue measure with density $\psi_{i}$ and $\Psi_{i}^{s}$ is singular to the Lebesgue measure with support $B\subset\mathbb{R}$.

Define the value function

[TABLE]

where $d\boldsymbol{\Psi}_{t}=\boldsymbol{\uppi}_{t-}d\boldsymbol{Q}_{t}-d\boldsymbol{\uppi}_{t}$, in which $\boldsymbol{\uppi}_{s-}=\boldsymbol{\uppi}$, $\boldsymbol{\uppi}_{t}=\boldsymbol{\uppi}\cdot\prodi_{[s,t]}(\boldsymbol{I}+d\boldsymbol{Q})(\boldsymbol{I}-d\boldsymbol{\Lambda})$ for any $t\in[s,T]$.

By the dynamic programming principle (see, for example, Fleming and Soner (2006)), for all $s\in[0,T]$ and $\Delta s>0$,

[TABLE]

which leads to

[TABLE]

If $s\notin B$, $\int_{\{s\}}d\Psi_{i}^{s}(u)=0$ for all $i\in\mathcal{X}$, and the Radon-Nikodym derivative $\psi_{i}$ exists at $s$, in which case the value function is a solution to the Hamilton-Jacobi-Bellman equations

[TABLE]

where

[TABLE]

as $\boldsymbol{\uppi}_{t-}(\boldsymbol{I}+d\boldsymbol{Q}_{t})d\boldsymbol{\Lambda}_{t}=d\boldsymbol{\Psi}_{t}$ and $\boldsymbol{\uppi}_{t-}=\boldsymbol{\uppi}\cdot\prodi_{[s,t)}(\boldsymbol{I}+d\boldsymbol{Q})(\boldsymbol{I}-d\boldsymbol{\Lambda})$, and

[TABLE]

If $s$ is an atomic point of $\Psi_{i}$, then $\pi_{i}(s-)>0$, and letting $\Delta s\rightarrow 0$ leads to the Hamilton-Jacobi-Bellman equations

[TABLE]

where $\boldsymbol{\Lambda}(\{s\})$ satisfies

[TABLE]

We claim that if a function $V$ satisfies (36) and (37), depending on whether $\boldsymbol{\Psi}$ is absolutely continuous or singular to the Lebesgue measure at $s$, then $V$ is the value function as defined in (35) for all waiting time $s\in[0,T]$ and distribution $\boldsymbol{\uppi}=(\pi_{1},\dots,\pi_{n})$ of initial states. To show this, we first prove $V$ is an upper bound of the value function. Since

[TABLE]

with (36) and (37), the left side of the above equation is no greater than

[TABLE]

with any $\boldsymbol{\Psi}$. Letting $s+\Delta s=T$ leads to

[TABLE]

which verifies that $V$ is an upper bound of the value function.

Now we show the equality can be achieved in (38) by providing the solution form of $\boldsymbol{\Psi}$ leading to $V$ that satisfies (36) and (37), and conclude that $V$ is indeed the value function.

From the Hamilton-Jacobi-Bellman equations (36),

[TABLE]

and the supremum on the right side is obtained by taking

[TABLE]

Note on the right side of (37)

[TABLE]

is a linear functional of $\boldsymbol{\uppi}\boldsymbol{\Lambda}(\{s\})$ with coefficient

[TABLE]

indicating the optimal $\boldsymbol{\Lambda}^{c}$ at any $s\in[0,T]$ should satisfy

[TABLE]

Combining the absolutely continuous and singular solutions (39) and (40) gives the desired solution form.

Finally, we show ${\eta}_{i}^{c}(s+)$ is indeed the expected life gain after waiting time $s$ if the occupancy is concentrated in state $i$ at $s-$ and none is transplanted at $s$. Let $\boldsymbol{e}_{i}$ denote the $n$-dimensional vector with the $i$th element being $1$ and all other elements being [math], then

[TABLE]

Since for any $t\in(s,T]$, $\boldsymbol{e}_{i}\cdot\prodi_{(s,t)}(\boldsymbol{I}+d\boldsymbol{Q})(\boldsymbol{I}-d\boldsymbol{\Lambda}^{c})$ is the limiting average occupancy at waiting time $t-$ if started from full occupancy of state $i$ at $s-$ and if none is transplanted at $s$, we have the desired result. The result for ${\gamma}_{i}^{c}(s+)$ can be shown similarly. ∎

Proof of Theorem 4.3

Proof.

By (18), it suffices to show that for all $s\in[0,T]$ and $i\in\mathcal{X}$,

[TABLE]

is a sufficient condition for

[TABLE]

where $\{\cdot\}_{i}$ indicates the $i$th element of the vector. Note for $j=1,2$,

[TABLE]

where $\{\boldsymbol{\Lambda}^{c_{j}}(t):t\in(s,T]\}$ are characterized in (16) and (17).

A key observation is that for any $s\in[0,T]$, the form of $\{\boldsymbol{\Lambda}^{c_{j}}(t):t\in(s,T]\}$ is independent of the initial condition $\boldsymbol{\uppi}_{s}$. Therefore for any $i\in\mathcal{X}$ and $s\in[0,T]$, $\{\boldsymbol{\Lambda}^{c_{j}}(t):t\in(s,T]\}$ also maximizes

[TABLE]

which is the objective under penalty parameter $c_{j}$ if only waiting time greater than $s$ is considered and all of the occupancy is concentrated at state $i$ at time $s$, and thus

[TABLE]

Subtracting $\{\boldsymbol{\upeta}^{c_{2}}(s+)-c_{2}\boldsymbol{\upgamma}^{c_{2}}(s+)\}_{i}$ from both sides of (41) gives

[TABLE]

The term on the right side of (42) is no less than $c_{2}-c_{1}$ as the probability of future transplantation

[TABLE]

Therefore

[TABLE]

and the desired result follows. ∎

Proof of Theorem 5.1

Proof.

The proof is similar to that of Theorem 3.3. For the multi-organ-type objective (28) define the value function

[TABLE]

The corresponding Hamilton-Jacobi-Bellman equations for the absolutely continuous case is

[TABLE]

where

[TABLE]

and

[TABLE]

and for the singular case is

[TABLE]

The forms of the optimal $\{\boldsymbol{\psi}^{(w)}_{s}:w=1,\dots,N_{o}\}$ and $\{\boldsymbol{\Psi}^{(w)}(\{s\}):w=1,\dots,N_{o}\}$ can be readily derived from the HJB equations and the constraints (26) and (27). ∎

Proof of Propositions 5.2 and 5.3

It suffices to show that if

[TABLE]

then

[TABLE]

where $c^{(w)}_{1}-c^{(w)}_{2}=\text{const}>0$ for all $w$ (Propositions 5.2), or $c^{(w)}_{1}>c^{(w)}_{2}$ and $c^{(v)}_{1}=c^{(v)}_{2}$ for all $v\neq w$ (Propositions 5.3).

It is obvious that (46) implies (48) with either condition on $\{c^{(w)}\}$. Similar to the proof of Theorem 4.3, we have

[TABLE]

Subtracting \sum_{w=1}^{N_{o}}\big{[}\eta_{i}^{(w),c_{2}}(s+)-c^{(w)}_{2}\gamma_{i}^{(w),c_{2}}(s+)\big{]} from both sides of (49) leads to the desired result, as the right side of the inequality becomes

[TABLE]