The Perfect Match: Pairing Organ Donors and Recipients using Multi-Objective Linear Optimization

Published in

Analytics Vidhya

6 min readJan 3, 2021

Image courtesy: https://www.fairtransplant.org/thelogo/

If you’ve ever been enthralled by the idea of having life after death, there’s a high probability of you encountering a person who is living a second life. Thanks to the phenomenal efforts made by medical researchers around the globe, that have given rise to advanced surgical procedures and finally made organ donation and transplantation a reality. However, despite having 90 percent of US adults support for organ donation, less than 60 percent of them actually sign up as a registered donor. Further statistics indicate the increasing rate of waiting list of recipients by 1 person every 9 minutes, with an average of 17 deaths per day the lives which could have been saved had they received the transplant at the right moment.

This article is an attempt to utilize the potential of large scale multi-objective optimization to devise a mechanism that efficiently pairs potential donors with their recipients ensuring maximum compatibility in a minimum timeframe and a high post-surgical survival rate.

All codes and data can be found via my GitHub repo: https://github.com/mohiteprathamesh1996/Heart-Donor-Pairing

Importing necessary dependencies

import random
import pandas as pd
from tqdm import tqdm
import matplotlib.pyplot as plt
import numpy as np
import itertools
import re
from pulp import *from IPython.core.display import display, HTMLdef display_side_by_side(dfs:list, captions:list):
    """Display tables side by side to save vertical space
    Input:
        dfs: list of pandas.DataFrame
        captions: list of table captions
    """
    output = ""
    combined = dict(zip(captions, dfs))
    for caption, df in combined.items():
        output += df.style.set_table_attributes("style='display:inline'").set_caption(caption)._repr_html_()
        output += "\xa0\xa0\xa0"
    display(HTML(output))import warnings
warnings.filterwarnings("ignore")

Data Description

For the purpose of this problem statement, I have created a sample dataset for recipients awaiting a heart transplant.

Pre-Surgical donor compatibility with recipient:

compatibility = pd.read_excel("donor_recepient_information.xlsx", sheet_name="pre_surgical_compatibility")compatibility.set_index(["Donor"], inplace=True)

This overall compatibility of heart transplant from Donor ‘i’ to Recipient ‘j’ is based on how well the donor’s heart fits comfortably inside the rib cage of the recipient.

2. Time to reach from donor to recipient (in hours):

time_to_reach = pd.read_excel("donor_recepient_information.xlsx", sheet_name="TimeToReach")time_to_reach.set_index(["Donor"], inplace=True)

The entire transplantation procedure must occur within 4 to 6 hours after a registered donor has officially died or has been declared brain dead. The above data indicates how much time (in hours) would it take to reach a prospective donor from it’s recipient.

3. Post-Surgical Survival Rate:

survival = pd.read_excel("donor_recepient_information.xlsx", sheet_name="post_surgical_survival_rate")survival.set_index(["Donor"], inplace=True)

This dataset reveals information about the projected survival rate of the recipient 4 to 5 years after receiving the transplant. Modern transplantation procedures have enabled to boost this metric to a range of 60–90%.

Decision Variables

Here we can define the binary decision variable, X(i, j) such that

X (i, j) = 1 if Donor ‘i’ is paired with Recipient ‘j’
and 0 otherwise, For all i ∈ [1,…,20] and j ∈ [1,…,20]

# List of donors
donor_list = compatibility.index.to_list()# List of recipients
recipient_list = compatibility.columns.to_list()# Dictionary of binary decision variables
var_dict = LpVariable.dicts(
name = "Match Found",
indexs = [(d, r) for d in donor_list for r in recipient_list], 
cat = "Binary")

Defining the Objective Functions

As mentioned earlier, our task now is to optimize 3 objective functions all at once i.e.:

Maximize -> Overall pre-surgical compatibility:

Maximize, Objective_1 = Σ Σ Compatibility(i, j) * X (i, j)
For all i ∈ [1,…,20] and j ∈ [1,…,20]

2. Minimize -> Overall time to reach potential donor

Minimize, Objective_2 = Σ Σ Time(i, j) * X (i, j)
For all i ∈ [1,…,20] and j ∈ [1,…,20]

3. Maximize -> Overall survival rate of all recipients post-surgery

Maximize, Objective_3 = Σ Σ Survival Rate(i, j) * X (i, j)
For all i ∈ [1,…,20] and j ∈ [1,…,20]

# Objective 1: Maximize donor compatibility
objective_1 = lpSum(
[compatibility.loc[d, r]*var_dict[(d, r)] \
for d in donor_list for r in recipient_list])
# Objective 2: Minimize time to reach the nearest donor
objective_2 = lpSum(
[time_to_reach.loc[d, r]*var_dict[(d, r)] \
for d in donor_list for r in recipient_list])
# Objective 3: Maximize post surgery survival probability
objective_3 = lpSum(
[survival.loc[d, r]*var_dict[(d, r)]\
 for d in donor_list for r in recipient_list])

When it comes to a multi-objective optimization problem, the most general technique would be to adopt a weighted sum approach, i.e. we basically assign weights to each of the objective functions and formulate a single-objective LP such as,

Maximize, Z = W1*Objective_1 + W2*(-Objective_2) + W3*Objective_3

Where W1, W2 and W3 are the weights assigned to each objective such that

W1, W2, W3 ∈ [0, 1] and W1 +W2 + W3

*In order to define the single objective function as a maximization problem, it’s important to have all independent objective functions in the same direction which is why you can see why I have multiplied the Objective_2, which is supposed to be a minimization problem, with a minus sign.

Constraints

The pairing must be done in a way such that each recipient can receive transplantation from exactly 1 donor and likewise every donor can donate to exactly 1 prospective recipient.

That is,

Σ X (i, j) = 1 where i (the donors) is a constant and j (the recipients)∈ [1,…,20]

Σ X (i, j) = 1 where j (the recipients) is a constant and i (the donors)∈ [1,…,20]

Obtaining the feasible space and corner point of optimality

Let us now formulate the above LP with the single objective function and including all constraints as a Python code and iterate over multiple combinations of weights and obtain a feasible space of optimal solutions. We will be running the code to obtain a new solutions for about 1000 iterations:

solutions = []for num in tqdm(range(1000)):
    weights = np.array([random.random(), random.random(), random.random()])# Ensuring sum of all weights adds up to 1
    weights = weights/sum(weights)# Define the maximization LPP
    model = LpProblem("Best Match", LpMaximize)# Define single objective
    model += (weights[0] * objective_1) + (weights[1] * -objective_2) + (weights[2] * objective_3)# LP constraints
    #----Constraint_1 -> A donor can donate to only one recipient
    for d in donor_list:
        model += lpSum([var_dict[(d, r)] for r in recipient_list]) == 1#----Constraint_2 -> A recipient can recieve from only 1 donor
    for r in recipient_list:
        model += lpSum([var_dict[(d, r)] for d in donor_list]) == 1# Solving the LP
    model.solve()# Saving results
    solutions.append([objective_1.value(), 
                      -objective_2.value(),
                      objective_3.value(), 
                      value(model.objective),
                      [(v.name,v.varValue) for v in model.variables() if v.varValue!=0]])# Save optimized results as a dataframe
df_optimal = pd.DataFrame(solutions, columns=["Objective_1", 
                                              "Objective_2",
                                              "Objective_3",
                                              "Single_Objective",
                                              "Optimal_Solution"])df_optimal.sort_values(by=["Single_Objective"], ascending=False).head()

Visualizing the 3-Dimensional feasible space

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt%matplotlib notebook
plt.rcParams["figure.figsize"] = (8, 5)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')ax.scatter(df_optimal[df_optimal["Single_Objective"]!=max(df_optimal["Single_Objective"])]["Objective_1"], 
           df_optimal[df_optimal["Single_Objective"]!=max(df_optimal["Single_Objective"])]["Objective_2"],
           df_optimal[df_optimal["Single_Objective"]!=max(df_optimal["Single_Objective"])]["Objective_3"], 
           c='green',
           marker='+',
           label = "Feasible solutions")ax.scatter(df_optimal[df_optimal["Single_Objective"]==max(df_optimal["Single_Objective"])]["Objective_1"], 
           df_optimal[df_optimal["Single_Objective"]==max(df_optimal["Single_Objective"])]["Objective_2"],
           df_optimal[df_optimal["Single_Objective"]==max(df_optimal["Single_Objective"])]["Objective_3"], 
           c='red',
           s=100,
           marker='*', 
           label = "Net Optimal Solution")ax.set_xlabel('Objective 1')
ax.set_ylabel('Objective 2')
ax.set_zlabel('Objective 3')plt.legend()plt.show()

We have finally arrived with our set of feasible solutions. The corner point (red star) is where you can find your optimal solution! The 3-dimensional surface formed as a result of all feasible solutions is also known as the Pareto Optimal Frontier named after the Italian engineer and economist Vilfredo Pareto (1848–1923).

Finding the best match

The ideal pair of the most compatible Donor s and Recipients can be obtained based on the optimal corner point in the above figure.

best_pair = df_optimal[df_optimal["Single_Objective"]==max(df_optimal["Single_Objective"])]status = pd.DataFrame(list(best_pair["Optimal_Solution"])[0], columns = ["Decisions","Compatibility"])status