{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Problem 13.1: Programming cellular response\n",
"\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"import eqtk"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have seen in this chapter that the combinatoric binding of multiple ligand types with multiple receptor types can lead to cellular addressing. In this problem, you will develop a scheme using two types of ligands, two types of \"A\" receptors, and two types of \"B\" receptors whereby the cell can respond to variation in ligand concentration in three orthogonal prescribed ways. These responses vary based on concentration of the respective receptors; the expression of these receptors delineate cell types. \n",
"\n",
"In our approach, we will closely follow the work of [Su, et al., 2022](https://doi.org/10.1016/j.cels.2022.03.001). We use the one-step promiscuous ligand-receptor model from this chapter, using the same notation. We define by a ligand \"word\" to be a set of specific concentration values for each ligand variant in a combination. For example, a concentration of 10 mM for ligand 1 and 100 mM for ligand 2 constitutes a ligand word (10 mM, 100 mM). In the diagram below, we have three ligand words denoted by the red circles, (1, 1000), (1000, 1), and (1000, 1000).\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"
\n",
"\n",
"Each word should encode a separate cellular response, giving different cell types. An example of three responses is shown in the image below.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"
\n",
"\n",
"The cell \"cares\" about the levels in the colored regions; the response in the hatched regions is irrelevant. In the above diagram, purple indicates low response and yellow indicates a high response.\n",
"\n",
"Your task is to come up with a set of parameters, eight $K_{ijk}$'s and eight $\\epsilon_{ijk}$'s, along with three sets of receptor concentrations, each containing the concentrations of A₁, A₂, B₁, and B₂, that have the desired response above, one each for each of the three words (1, 1000), (1000, 1), and (1000, 1000). An example of a desired response is shown below.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"
\n",
"\n",
"In working this problem, it helps to have some handy functions from the chapter."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def make_rxns(nA, nB, nL):\n",
" '''\n",
" Generate trimolecular binding reactions for a system with\n",
" nA, nB, and nL types of Type A receptors, Type B receptors,\n",
" and ligands, respectively. Returns a single string.\n",
" '''\n",
" rxns = \"\"\n",
" for k in range(nB):\n",
" for j in range(nL):\n",
" for i in range(nA):\n",
" rxns += f\"A_{i+1} + L_{j+1} + B_{k+1} <=> T_{i+1}_{j+1}_{k+1}\\n\"\n",
" return rxns\n",
"\n",
"\n",
"def make_N(nA, nB, nL):\n",
" \"\"\"Generate a stoichiometric matrix for ligand-receptor binding with\n",
" nA, nB, and nL types of Type A receptors, Type B receptors, and \n",
" ligands, respectively.\n",
" \"\"\"\n",
" rxns = make_rxns(nA, nB, nL)\n",
" N = eqtk.parse_rxns(rxns)\n",
"\n",
" # Sorted names\n",
" names = sorted(N.columns, key=lambda s: (len(s), s))\n",
"\n",
" # Sorted columns\n",
" N = N[names]\n",
"\n",
" # As a Numpy array\n",
" return N.to_numpy(copy=True, dtype=float)\n",
"\n",
"\n",
"def readout(epsilon, c):\n",
" \"\"\"Readout function of a set of complex concentrations. This\n",
" is the expression level of the genes under regulation of\n",
" singaling from the receptors.\n",
" \"\"\"\n",
" return np.dot(epsilon, c[:, -len(epsilon):].transpose())"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A_1 + L_1 + B_1 <=> T_1_1_1\n",
"A_2 + L_1 + B_1 <=> T_2_1_1\n",
"A_1 + L_2 + B_1 <=> T_1_2_1\n",
"A_2 + L_2 + B_1 <=> T_2_2_1\n",
"A_1 + L_1 + B_2 <=> T_1_1_2\n",
"A_2 + L_1 + B_2 <=> T_2_1_2\n",
"A_1 + L_2 + B_2 <=> T_1_2_2\n",
"A_2 + L_2 + B_2 <=> T_2_2_2\n",
"\n"
]
}
],
"source": [
"print(make_rxns(2, 2, 2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can make a stoichiometric matrix using these functions; it will be the same for all calculations we do."
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"N = make_N(2, 2, 2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Half of the battle in problems like these is keeping track of indices to know what species is what. The columns of the stoichiometric matrix correspond to the following species.\n",
"\n",
"| index | species |\n",
"| ----------- | ----------- |\n",
"| 0 | A₁ |\n",
"| 1 | A₂ |\n",
"| 2 | B₁ |\n",
"| 3 | B₂ |\n",
"| 4 | L₁ |\n",
"| 5 | L₂ |\n",
"| 6 | A₁L₁B₁ |\n",
"| 7 | A₂L₁B₁ |\n",
"| 8 | A₁L₂B₁ |\n",
"| 9 | A₂L₂B₁ |\n",
"| 10 | A₁L₁B₂ |\n",
"| 11 | A₂L₁B₂ |\n",
"| 12 | A₁L₂B₂ |\n",
"| 13 | A₂L₂B₂ |\n",
"\n",
"These indices also correspond to the species in our inputted initial concentrations and in the output of EQTK's solver for the equilibrium concentrations. With that in mind, we can write a function to generate our initial concentration values for a given set of initial ligand concentrations, `cL0`."
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"def make_c0_grid(cL0, cA10, cA20, cB10, cB20):\n",
" \"\"\"Create an array of initial concentrations.\n",
" \n",
" Parameters\n",
" ----------\n",
" cL0 : array_like\n",
" Array of ligand concentrations. This is the same for ligand\n",
" 1 and ligand 2, so a single 1D array is all that is required.\n",
" cA10 : float\n",
" Total concentration of receptor A1.\n",
" cA20 : float\n",
" Total concentration of receptor A2.\n",
" cB10 : float\n",
" Total concentration of receptor A1.\n",
" cB20 : float\n",
" Total concentration of receptor B2.\n",
"\n",
" Returns\n",
" -------\n",
" output : 2D Numpy array\n",
" If n is the length of the cL0 input array, the output\n",
" is an n² by 14 array. Columns 0 through 3 are the total\n",
" receptor concentrations. Columns 4 and 5 are the ligand\n",
" concentrations. Columns 6 through 13 are the concentrations\n",
" of all of the possible trimers (all set to zero).\n",
" \"\"\"\n",
" # Number of ligand concentrations\n",
" n = len(cL0)\n",
" \n",
" # Ligand concentrations\n",
" cL0 = np.meshgrid(*tuple([cL0] * 2))\n",
"\n",
" # Initialize c0\n",
" c0 = np.zeros((n**2, 14))\n",
"\n",
" # Add ligand concentrations\n",
" c0[:, 4] = cL0[0].flatten()\n",
" c0[:, 5] = cL0[1].flatten()\n",
" \n",
" # Add receptor concentrations\n",
" c0[:, 0] = cA10\n",
" c0[:, 1] = cA20\n",
" c0[:, 2] = cB10\n",
" c0[:, 3] = cB20\n",
"\n",
" return c0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In trying to find parameters, we will use the ligand concentrations below."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"cL0 = np.array([1.0, 32.0, 1000.0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, in order to do the calculation, we need to specify the prescribed target response of the various cell types and for which ligand words the target response is relevant (ignoring the hatches areas in the graphic above). We refer to the latter as `active_words`, defined below."
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"# Ligand words from top diagram\n",
"active_words = np.array([0, 0, 1, 0, 0, 0, 1, 0, 1], dtype=bool)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can see how it corresponds to the words by comparing to our target ligand concantrations."
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Ligand 1: [ 1. 32. 1000. 1. 32. 1000. 1. 32. 1000.]\n",
"Ligand 2: [ 1. 1. 1. 32. 32. 32. 1000. 1000. 1000.]\n",
"active words: [False False True False False False True False True]\n"
]
}
],
"source": [
"c0 = make_c0_grid(np.array([1, 32, 1000]), 0, 0, 0, 0)\n",
"print(\"Ligand 1: \", c0[:, 4])\n",
"print(\"Ligand 2: \", c0[:, 5])\n",
"print(\"active words:\", active_words)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now specify a list of target responses."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"targets = [\n",
" np.array([0, 0, 0, 0, 0, 0, 1, 0, 0]), # Ligand 1 low, ligand 2 high\n",
" np.array([0, 0, 1, 0, 0, 0, 0, 0, 0]), # Ligand 1 high, ligand 2 low\n",
" np.array([0, 0, 0, 0, 0, 0, 0, 0, 1]), # Ligand 1 high, ligand 2 high\n",
"]\n",
"\n",
"# Useful to know how many cell types we have\n",
"n_cell_types = len(targets)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, it will be useful to have a function to solve for the normalized response of a cell to an inputted set of initial concentrations and parameters. The function returns the cellular response, with the maximum value of the response over all ligand concentrations being set to one. That is, we want a function that returns $s_\\mathrm{norm}$, defined by\n",
"\n",
"\\begin{align}\n",
"&s(l_1, l_2) = \\sum_{ijk} \\epsilon_{ijk}\\,t_{ijk}(l_1, l_2),\\\\[1em]\n",
"&s_\\mathrm{norm}(l_1,l_2) = \\frac{s(l_1, l_2)}{\\text{maximal }s\\text{ over all ligand concentrations}}.\n",
"\\end{align}\n",
"\n",
"**a)** Complete the function below."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def solve_norm_response(c0, N, K, epsilon):\n",
" \"\"\"Solve for normalized response. Returns a one-dimensional array\n",
" containing signal levels corresponding to each row in `c0`.\n",
" \"\"\"\n",
" pass"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**b)** To find optimal parameters that give cell responses close to the target responses, you can use whatever method you like. However, we will walk you through a method following [Su, et al.](https://doi.org/10.1016/j.cels.2022.03.001), which you may find easier to do that developing your own method.\n",
"\n",
"We will use a least squares method. To do this we define a **residual**, which is the value of the signal a given active ligand word, minus the target signal for that word. \n",
"\n",
"\\begin{align}\n",
"\\text{residual}(l_1, l_2) = s_\\mathrm{norm}(l_1, l_2) - s_\\mathrm{norm}^\\mathrm{target}(l_1, l_2).\n",
"\\end{align}\n",
"\n",
"We then seek to minimize the sum of these residuals. Note that there exists three residuals for each cell type in the above example, for a total of nine residuals. \n",
"\n",
"The `scipy.optimize.least_squares()` function performs this calculation automatically. Its call signature is\n",
"\n",
" scipy.optimize.least_squares(residuals, x0, bounds=(lower_bounds, upper_bounds), args=args)\n",
"\n",
"Here, `residuals()` is a function that returns the array of residuals and has call signature\n",
"\n",
" residuals(x, *args),\n",
"\n",
"where `x` is an array of the parameters being optimized and `args` is a tuple of any other arguments the function needs. To keep things organized, we will take the first eight entries in `x` to be the $K_{ijk}$ values, the next eight entries to be the $\\varepsilon_{ijk}$ values, and the remaining twelve entries to be the receptor concentrations for the respective cell types.\n",
"\n",
"Complete the function below. "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def residuals(x, cL0, targets, active_words, N):\n",
" \"\"\"Return array of residuals.\"\"\"\n",
" # Unpack x, log K's, then log epsilons, then receptor concs\n",
" K = x[:8]\n",
" epsilon = x[8:16]\n",
" receptor_conc = x[16:]\n",
"\n",
" pass"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given our function residuals, we can define `args` to match its call signature."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"args = (cL0, targets, active_words, N)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**c)** To complete the optimization, we need to two more pieces. First, we need to provide bounds for all of the parameters. All of the $K_{ijk}$, $\\varepsilon_{ijk}$, and receptor concentrations must be positive, so we have the following for the lower bounds."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"lower_bounds = np.zeros(16 + 4 * n_cell_types)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the upper bounds, Su, et al., specified that all $K_{ijk}$ and $\\varepsilon_{ijk}$ must be between zero and one. This effectively sets the arbitrary units of the ligand and receptor concentrations. (We do not require that the sum of the $K_{ijk}$'s and the sum of the $\\varepsilon_{ijk}$'s sum to one.) The receptor concentrations are unbounded, so we can set the upper bounds, recalling the ordering of our parameter array."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"upper_bounds = np.concatenate(\n",
" (\n",
" np.ones(16), # Upper bounds on K's and ε's\n",
" np.inf * np.ones(4 * n_cell_types) # Upper bounds on receptor concs\n",
" )\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we must provide a guess of the parameter values. Su, et al. used the strategy that the initial guess of all receptor concentrations were one and randomly chose values for $K_{ijk}$ and $\\epsilon_{ijk}$ from a uniform distribution on the interval \\[0, 1\\]. You can try this strategy as in the code cell below or adopt one of your own."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# Initial guesses of parameter values\n",
"K0 = np.random.uniform(size=8)\n",
"epsilon0 = np.random.uniform(size=8)\n",
"receptor_conc0 = np.ones(4 * n_cell_types)\n",
"\n",
"# Package together\n",
"x0 = np.concatenate((K0, epsilon0, receptor_conc0))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, you can perform the optimization be running\n",
"\n",
" res = scipy.optimize.least_squares(\n",
" residuals, x0, bounds=(lower_bounds, upper_bounds), args=args\n",
" )\n",
" \n",
"The object return is the optimization result. You can pull out the parameter values from the `res.x` attribute.\n",
"\n",
" K = res.x[:8]\n",
" epsilon = res.x[8:16]\n",
" receptor_conc = res.x[16:]\n",
" \n",
"Find a set of parameters that give the desired cell types. Make a plot like the triptych above for the various cell types to demonstrate that you parameter set does indeed give the desired orthogonal cell types. Also report the respective receptor concentration profiles of each of the three cell types.\n",
"\n",
"*Hint:* The optimization will converge to a local minimum, which may not hit your target addressing very well. There are *lots* of local minima. You can try rerunning the optimization with a different set of starting guesses for the parameters and receptor concentrations. You may have to do this several times."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**d)** If you like, try coming up with parameter values for other word sets and/or target responses."
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [],
"source": [
"c0\n",
"\n",
"# Set up fixed concentrations\n",
"fixed_c = -np.ones_like(c0)\n",
"fixed_c[:, 4:6] = c0[:, 4:6]\n",
"\n",
"# # Solve for the concentrations\n",
"# c = eqtk.fixed_value_solve(c0=c0, fixed_c=fixed_c, N=N, K=K)\n",
"\n",
"# # Build and return normalized response\n",
"# s = readout(epsilon, c)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(9, 14)"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"c0.shape"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ -1., -1., -1., -1., 1., 1., -1., -1., -1.,\n",
" -1., -1., -1., -1., -1.],\n",
" [ -1., -1., -1., -1., 32., 1., -1., -1., -1.,\n",
" -1., -1., -1., -1., -1.],\n",
" [ -1., -1., -1., -1., 1000., 1., -1., -1., -1.,\n",
" -1., -1., -1., -1., -1.],\n",
" [ -1., -1., -1., -1., 1., 32., -1., -1., -1.,\n",
" -1., -1., -1., -1., -1.],\n",
" [ -1., -1., -1., -1., 32., 32., -1., -1., -1.,\n",
" -1., -1., -1., -1., -1.],\n",
" [ -1., -1., -1., -1., 1000., 32., -1., -1., -1.,\n",
" -1., -1., -1., -1., -1.],\n",
" [ -1., -1., -1., -1., 1., 1000., -1., -1., -1.,\n",
" -1., -1., -1., -1., -1.],\n",
" [ -1., -1., -1., -1., 32., 1000., -1., -1., -1.,\n",
" -1., -1., -1., -1., -1.],\n",
" [ -1., -1., -1., -1., 1000., 1000., -1., -1., -1.,\n",
" -1., -1., -1., -1., -1.]])"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fixed_c"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.10"
}
},
"nbformat": 4,
"nbformat_minor": 4
}