Python Pulp Integer Linear Program with dynamic constraint

I want to solve a mixed integer linear program with the following objective function:

J = maximize (f1(x) + f2(x)) subject to constraint: cost(x) <= threshold

where x is the set of selected variables, f1 and f2 are two scoring functions and cost is the cost function.

f2 is a function based on similarity between the selected variables. I don't know how to formulate this function in pulp.

This is my minimal working example in which function f2 is the similarity between two ingredients and I want to add similarity[i][j] to the objective function if j is already in selected variables, but don't know how to do it.

import numpy as np
import pulp
threshold = 200
model = pulp.LpProblem('selection', pulp.LpMaximize)
similarity = np.array([[1., 0.08333333, 0.1, 0., 0., 0.0625],
                       [0.08333333, 1., 0.33333333,
                           0., 0.11111111, 0.07692308],
                       [0.1, 0.33333333, 1., 0.2, 0., 0.09090909],
                       [0., 0., 0.2, 1., 0., 0.],
                       [0., 0.11111111, 0., 0., 1., 0.27272727],
                       [0.0625, 0.07692308, 0.09090909, 0., 0.27272727, 1.]])
ingredients = ['var_%d' % i for i in range(6)]
scores = np.random.randint(1, 3, size=len(ingredients))
costs = np.random.randint(20, 60, len(ingredients))
scores = dict(zip(ingredients, scores))
costs = dict(zip(ingredients, costs))
x = pulp.LpVariable.dict(
    'x_%s', ingredients, lowBound=0, upBound=1, cat=pulp.LpInteger)
model += sum([scores[i] * x[i] for i in ingredients])
model += sum([costs[i] * x[i] for i in ingredients]) <= threshold
solver = pulp.solvers.PULP_CBC_CMD()
model.solve(solver)

This code basically considers only static costs (encoded in costs variable). How can I dynamically add similarity costs that are the similarity variable?