4D Heatmap for PDP

Right now I'm working on a custom implementation of a PDP interaction for 3 features at the same time and I have a problem with visualizing the data. Sometime ago I thought that it would be nice to represent a 3D-heatmap for features' interaction, where the surface shape represents PDP interaction for 2 out of 3 features (x&y and x&z), and the heatmap texture layer represents interaction of the last two features (z&y).

Right now I have a working (barebones, though) implementation of the PDP heatmap for 2 features.

def pdp_custom_3D(estimated_model, X, y, n_splits, target_name_1, target_name_2, prefit = True, X_train = None, y_train = None):

import matplotlib.pyplot as plt
import seaborn as sb

if prefit == False:
    try:
        estimated_model.fit(X_train, y_train)
    except:
        logging.warning("Estimated model must have fit method.")
            
PDP_list_1 = list()
PDP_list_2 = list()

feature_list_1 = list()
feature_list_2 = list()

x_max_1 = X[target_name_1].max()
x_min_1 = X[target_name_1].min()

x_max_2 = X[target_name_2].max()
x_min_2 = X[target_name_2].min()


print('x_max_1', x_max_1)
print('x_min_1', x_min_1)
print('x_max_2', x_max_2)
print('x_min_2', x_min_2)

X_copy = X.copy() 

step_1 = abs(x_max_1 - x_min_1)/n_splits
step_2 = abs(x_max_2 - x_min_2)/n_splits

x_axis = int(n_splits+1)
y_axis = int(n_splits+1)

s = (x_axis, y_axis)  
PDP_result = np.zeros(s)

counter_1 = x_min_1
counter_2 = x_min_2


for i in range(n_splits+1):
    
    feature_list_1.append(counter_1 + (i*step_1))
    X_copy[target_name_1] = counter_1 + i * step_1

    for j in range(n_splits+1):
        
        feature_list_2.append(counter_2 + j * step_2)
        X_copy[target_name_2] = counter_2 + j * step_2
        temp = estimated_model.predict(X_copy)
        PDP_result[i][j] = temp.mean()
    
sb.heatmap(PDP_result)


return  
    
    

And It gets me somewhat nice heatmap

Is there any way I can get 3D heatmap (3 axis are selected features and 4th is mean of the prediction model is giving me back - It is being represented as a heat signature here, and I want a surface as a second representation of the prediction) with the same data forms?

Somewhat like this

Thank you all and have a nice day!