Rat model: implementing with real data

Rationale

When trying to model the D'Aquila et al. (2012) data, I ran into an additional problem. How do we actually plug in this real-world data into our model in a way that makes sense?

Here I describe how I can't model the D'Aquila (2012) data in a way that makes sense with the existing model, and I'm proposing a modified solution.



In [1]:

    
from minimal_example_interface import *



In [2]:

    
def multiplicative_term(kappa_val,r_array,cue_array):
    assert(type(kappa_val)==float or type(kappa_val)==int)
    assert(type(r_array)==np.ndarray)
    assert(type(cue_array)==np.ndarray)
    
    return({"kappa": kappa_val, "r": r_array,"cue": cue_array})

set_term_function(multiplicative_term)



In [3]:

    
def SaltSugarModelWithCue(term_list):
    return sum([SaltSugarModelWithCue_SingleTerm(term) for term in term_list])

set_model_calculation_function(SaltSugarModelWithCue)



In [4]:

    
def SaltSugarModelWithCue_SingleTerm(term):
    return term['kappa']*term['r']*term['cue']

set_singleterm_calculation_function(SaltSugarModelWithCue_SingleTerm)

Background

The model we proposed in our draft was

\begin{align} \begin{split} \tilde r(\mathbf{r},\mathbf{\kappa})= & \kappa_{\text{Na}} r_{\text{Na}}\mathbf{c}+\kappa_{h} r_{h}\mathbf{c} \end{split} \end{align}

where

the two terms $\mathbf{\kappa}_{\text{Na}}$ and $\mathbf{\kappa}_{h}$ are the "salt need" and "hypernatremia avoidance need".
the two terms $\mathbf{r}_{\text{Na}}$ and $\mathbf{r}_{h}$ are the "salt expected value" and "hypernatremia expected value" for the salt solution.
$\mathbf{c}$ is the cue value which we can ignore in this context because it is constant
$\tilde r(\mathbf{r},\mathbf{\kappa})$ is the incentive salience of the solution.

Under salt deprivation these are:

\begin{align} \tilde r(\mathbf{r},\mathbf{\kappa})= & \mathbf{\kappa}_{\textrm{Na}} \mathbf{r}_\textrm{Na}+\mathbf{\kappa}_{h} \mathbf{r}_{h} \\ \nonumber = & 1 \times \begin{bmatrix}0.5 \\ 1.0\end{bmatrix} + 2 \times \begin{bmatrix}0.0 \\ -1.0\end{bmatrix} \\ \nonumber = & \begin{bmatrix}0.5 \\ -1.0\end{bmatrix} \begin{matrix}\textit{Moderate Salt solution} \\ \textit{Strong Salt solution}\end{matrix} \end{align}

To model salt satiation, we all we have to do is flip the sign on Salt Need:

\begin{align} \tilde r(\mathbf{r},\mathbf{\kappa})= & \mathbf{\kappa}_{\textrm{Na}} \mathbf{r}_\textrm{Na}+\mathbf{\kappa}_{h} \mathbf{r}_{h} \\ \nonumber = & -1 \times \begin{bmatrix}0.5 \\ 1.0\end{bmatrix} + 2 \times \begin{bmatrix}0.0 \\ -1.0\end{bmatrix} \\ \nonumber = & \begin{bmatrix}-0.5 \\ -3.0\end{bmatrix} \begin{matrix} \textit{Moderate Salt solution}\\ \textit{Strong Salt solution} \end{matrix} \end{align}

These were formed out of the following objectives:

Under at least weak to moderate salt depletion, the rat prefers the moderate solution to the strong solution; it may approach the strong solution sometimes, but less often than the moderate solution. So under salt deprivation, the solutions will be actively approached, but the strong solution will be less preferable than the moderate solution.
Under salt satiation, the rat prefers the moderate solution to the strong solution (i.e., finds it less aversive). So the preference ordering from (1) is maintained.
Under sufficiently strong salt depletion, the strong solution can become appetitive while still being less appetitive than the moderate solution.
The rat's choice of the moderate or salty solution will be a function of the relative incentive salience or "wanting" of the two solutions, where incentive salience is a function of cue strength for the two solutions. Thus, under salt depletion, the rat will take the strong solution if there is no cue for the moderate solution.

So how do we incorporate the real-world data we have in D'Aquilla (2012)?

Incorporating the real-world data

The dataset contains:

Strength of the salt solution, which was varied from one condition to the next
Salt deprivation
A variety of outcome measures, including number of licks and latency until first lick of the stimulus.

So, considering our terms:

\begin{align} \begin{split} \tilde r(\mathbf{r},\mathbf{\kappa})= & \kappa_{\text{Na}} r_{\text{Na}}+\kappa_{h} r_{h} \end{split} \end{align}

how should the data be represented in the equation?

Well, first, we can represent the "latency until first lick" as $\tilde r(\mathbf{r},\mathbf{\kappa})$ - so we're trying to predict that outcome from the numbers.

We could represent strength of the salt solution as the expected reward of the salt ($r_{\text{Na}}$). We can do this because the rats have already learned to associate stimuli with a certain level of salt.

The D'Aquilla (2012) data is actually very straightforward. You can represent it with only a single term. The salt levels are just "low" (0.09%) and "moderate" (2.7%), not "high" like the dead-sea-salt model (8.9%).

Incorporating additional data

So I want to incorporate hypothetical data, as a stand-in until I can get some empirical data on 8.9%.

So now we're trying to model three solutions, with values 0.009, 0.027, and 0.089, using the term $r_{\text{Na}}$.

But if we posit that the expected value for satiating salt need is the stimulus strength itself, we should also represent the hypertremia danger $r_{h}$ using the same number. That's not going to work in a straightforward way, because in the model as previously defined, I used $[1.0, 1.0]$ to represent the $r_{\text{Na}}$ for the two stimuli, but $[0.1, 1.0]$ to represent the $r_{h}$ for the two stimuli. Now we're using the exact same values.

$r_{\text{h}}=r_{\text{Na}}$ is going to be a problem for us. But we could represent the non-linear dynamics above with a simple: $r_{\text{h}}=(r_{\text{Na}})^2$.

So writing those values out:

$r_{\text{Na}} = [0.009, 0.027, 0.89]$

$r_{\text{h}} = (r_{\text{Na}})^2 = [0.000081, 0.000729, 0.007921]$

We can still modify salt need and hypernatremia parameters arbitrarily to get a best fit. So this turns out to give us $\tilde r(\mathbf{r},\mathbf{\kappa})$ values that check out.

But I worry people are going to ask about a justification for $r_{\text{h}}=(r_{\text{Na}})^2$. My justification is that hypernatremia danger is very low for low to moderate levels of salt, but increases exponentially as you are drinking higher and higher levels of salt concentration. Intuitively that sounds right, but I don't know how to justify it any better than simply asserting it.

The code below implements the model in python.

We can then create variables to represent each term and then calculate their value:



In [9]:

    
cue=np.array([1.0,1.0,1.0])
term_Na=multiplicative_term(10,np.array([0.009, 0.027, 0.089]),cue)
term_h=multiplicative_term(-100,np.array([pow(x,2) for x in [0.009, 0.027, 0.089]]),cue)
term_Glc=multiplicative_term(1,np.array([0.0,0.0,0.0]),cue)

We can write this out mathematically, first under salt depletion:



In [10]:

    
Markdown(get_and_format_main_equation_text(term_Na,term_h,term_Glc,three_solutions=['Weak Salt solution','Moderate salt solution','Strong salt solution']))









    Out[10]:




$$
\begin{align}
\begin{split}
\tilde r(\mathbf{r},\mathbf{\kappa},\mathbf{c})= & 
\kappa_{\text{Na}} \mathbf{r}_{\text{Na}}\mathbf{c} 
& + \kappa_{h} \mathbf{r}_{h}\mathbf{c}
& + \kappa_{\text{Glc}} \mathbf{r}_{\text{Glc}}\mathbf{c} \\
 = & 10 \times \begin{bmatrix}0.01 \\ 0.03 \\ 0.09\end{bmatrix}\circ \begin{bmatrix}1.0 \\ 1.0 \\ 1.0\end{bmatrix} & + -100 \times \begin{bmatrix}0.0001 \\ 0.0007 \\ 0.0079\end{bmatrix}\circ \begin{bmatrix}1.0 \\ 1.0 \\ 1.0\end{bmatrix} & + 1 \times \begin{bmatrix}0.0 \\ 0.0 \\ 0.0\end{bmatrix}\circ \begin{bmatrix}1.0 \\ 1.0 \\ 1.0\end{bmatrix} \\
 = & \begin{bmatrix}0.09 \\ 0.27 \\ 0.89\end{bmatrix} & + \begin{bmatrix}-0.0081 \\ -0.0729 \\ -0.7921\end{bmatrix} & + \begin{bmatrix}0.0 \\ 0.0 \\ 0.0\end{bmatrix} \\
 = & \begin{bmatrix}0.08 \\ 0.2 \\ 0.1\end{bmatrix} \begin{matrix} \textit{Weak Salt solution}\\ \textit{Moderate salt solution}\\ \textit{Strong salt solution} \end{matrix}
\end{split}
\end{align}
$$

These are in the right order. The weak salt solution is least preferred (in accordance with the D'Aquilla, 2012 data), the moderate solution is next preferred, and the strong solution is not as preferred as the moderate solution.

Now, under salt satiation (just flipping the sign of the salt need variable)



In [8]:

    
term_Na=multiplicative_term(-10,np.array([0.009, 0.027, 0.089]),cue)

Markdown(get_and_format_main_equation_text(term_Na,term_h,term_Glc,three_solutions=['Weak Salt solution','Moderate salt solution','Strong salt solution']))









    Out[8]:




$$
\begin{align}
\begin{split}
\tilde r(\mathbf{r},\mathbf{\kappa},\mathbf{c})= & 
\kappa_{\text{Na}} \mathbf{r}_{\text{Na}}\mathbf{c} 
& + \kappa_{h} \mathbf{r}_{h}\mathbf{c}
& + \kappa_{\text{Glc}} \mathbf{r}_{\text{Glc}}\mathbf{c} \\
 = & -10 \times \begin{bmatrix}0.01 \\ 0.03 \\ 0.09\end{bmatrix}\circ \begin{bmatrix}1.0 \\ 1.0 \\ 1.0\end{bmatrix} & + -100 \times \begin{bmatrix}0.0001 \\ 0.0007 \\ 0.0079\end{bmatrix}\circ \begin{bmatrix}1.0 \\ 1.0 \\ 1.0\end{bmatrix} & + 1 \times \begin{bmatrix}0.0 \\ 0.0 \\ 0.0\end{bmatrix}\circ \begin{bmatrix}1.0 \\ 1.0 \\ 1.0\end{bmatrix} \\
 = & \begin{bmatrix}-0.09 \\ -0.27 \\ -0.89\end{bmatrix} & + \begin{bmatrix}-0.0081 \\ -0.0729 \\ -0.7921\end{bmatrix} & + \begin{bmatrix}0.0 \\ 0.0 \\ 0.0\end{bmatrix} \\
 = & \begin{bmatrix}-0.1 \\ -0.34 \\ -1.68\end{bmatrix} \begin{matrix} \textit{Weak Salt solution}\\ \textit{Moderate salt solution}\\ \textit{Strong salt solution} \end{matrix}
\end{split}
\end{align}
$$

Here, all solutions are negative (aversive), and the strong solution is the most aversive of them all.

If we move ahead with this, then I can run a formal model, fine-tuning the $\kappa$ values by subject to find a best fit for the precise behavior shown for each subject.

An alternative model

But I worry that the reviewers won't like my arbitrary non-linear transformation $r_{\text{h}}=(r_{\text{Na}})^2$. There is an alternative, but the model gets more complicated. This is to apply an inverse logit function:

Then we can actually use the same values for $r_{\text{h}}$ as for $r_{\text{Na}}^2$, for a high salt need:

\begin{align} \tilde r(\mathbf{r},\mathbf{\kappa})= & \text{mean}( \text{logit}^{-1}(\mathbf{\kappa}_{\textrm{Na}} \mathbf{r}_\textrm{Na}),\text{logit}^{-1}(\mathbf{\kappa}_{h} \mathbf{r}_{h}))) \\ \nonumber = & \text{mean}(\text{logit}^{-1}( 100 \times \begin{bmatrix}0.009 \\ 0.027 \\ 0.089 \end{bmatrix}) ,\text{logit}^{-1}( -50 \times \begin{bmatrix}0.009 \\ 0.027 \\ 0.089 \end{bmatrix} )) \\ \nonumber = & \text{mean}( \begin{bmatrix}0.71 \\ 0.94 \\ 1.00 \end{bmatrix} , \begin{bmatrix}0.39 \\ 0.21 \\ 0.01 \end{bmatrix} ) \\ \nonumber = & \begin{bmatrix}0.55 \\ 0.57 \\ 0.51 \end{bmatrix} \begin{matrix} \textit{Weak salt solution} \\ \textit{Moderate Salt solution} \\ \textit{Strong Salt solution} \end{matrix} \end{align}

...and for low salt need (just flipping the sign)

\begin{align} \tilde r(\mathbf{r},\mathbf{\kappa})= & \text{mean}( \text{logit}^{-1}(\mathbf{\kappa}_{\textrm{Na}} \mathbf{r}_\textrm{Na}),\text{logit}^{-1}(\mathbf{\kappa}_{h} \mathbf{r}_{h}))) \\ \nonumber = & \text{mean}(\text{logit}^{-1}( -100 \times \begin{bmatrix}0.009 \\ 0.027 \\ 0.089 \end{bmatrix}) ,\text{logit}^{-1}( -50 \times \begin{bmatrix}0.009 \\ 0.027 \\ 0.089 \end{bmatrix} )) \\ \nonumber = & \text{mean}( \begin{bmatrix}0.29 \\ 0.06 \\ 0.00 \end{bmatrix} , \begin{bmatrix}0.39 \\ 0.21 \\ 0.01 \end{bmatrix} ) \\ \nonumber = & \begin{bmatrix}0.34 \\ 0.13 \\ 0.01 \end{bmatrix} \begin{matrix} \textit{Weak salt solution} \\ \textit{Moderate Salt solution} \\ \textit{Strong Salt solution} \end{matrix} \end{align}

If you like, we can then apply a logit function to 'undo' the inverse logit function. Values higher than 0.5 will be positive and values lower than 0.5 will be negative:

\begin{align} \text{logit}(\begin{bmatrix}0.55 \\ 0.57 \\ 0.51 \end{bmatrix} ) = \begin{bmatrix}0.20 \\ 0.29 \\ 0.02 \end{bmatrix} \begin{matrix} \textit{Weak salt solution} \\ \textit{Moderate Salt solution} \\ \textit{Strong Salt solution} \end{matrix} \end{align}\begin{align} \text{logit}(\begin{bmatrix}0.34 \\ 0.13 \\ 0.01 \end{bmatrix} ) = \begin{bmatrix} -0.67 \\ -1.86 \\ -5.14 \end{bmatrix} \begin{matrix} \textit{Weak salt solution} \\ \textit{Moderate Salt solution} \\ \textit{Strong Salt solution} \end{matrix} \end{align}

Again, as with the model above, this accomplishes all our objectives above. At high salt need, the moderate solution is the most attractive, followed by the weak solution and then the strong solution. At low salt need, the sign flips to negative (when applying the inverse logit function), and the weak salt solution is the least aversive, followed by the moderate and then strong solution being most aversive.

Next steps

The "real world" data I have incorporated here is actually just the salt solutions! What we do have that we didn't have is D'Aquilla (2012) data confirming that weak and moderate solutions work the way Berridge described.

We can also get (from Berridge's team) data showing that the strong solution is more aversive under salt satiation. But that doesn't help us with what really need: showing relative preferences for moderate vs. strong solutions. So if we really want to talk about "real world data" we need that.

But I'm hoping to get your thoughts on these two model variants.