Systems biology reveals uncoupling beyond UCP1 in human white fat-derived beige adipocytes

Pharmaceutical induction of metabolically active beige adipocytes in the normally energy storing white adipose tissue has potential to reduce obesity. Mitochondrial uncoupling in beige adipocytes, as in brown adipocytes, has been reported to occur via the uncoupling protein 1 (UCP1). However, several previous in vitro characterizations of human beige adipocytes have only measured UCP1 mRNA fold increase, and assumed a direct correlation with metabolic activity. Here, we provide an example of pharmaceutical induction of beige adipocytes, where increased mRNA levels of UCP1 are not translated into increased protein levels, and perform a thorough analysis of this example. We incorporate mRNA and protein levels of UCP1, time-resolved mitochondrial characterizations, and numerous perturbations, and analyze all data with a new fit-for-purpose mathematical model. The systematic analysis challenges the seemingly obvious experimental conclusion, i.e., that UCP1 is not active in the induced cells, and shows that hypothesis testing with iterative modeling and experimental work is needed to sort out the role of UCP1. The analyses demonstrate, for the first time, that the uncoupling capability of human beige adipocytes can be obtained without UCP1 activity. This finding thus opens the door to a new direction in drug discovery that targets obesity and its associated comorbidities. Furthermore, the analysis advances our understanding of how to evaluate UCP1-independent thermogenesis in human beige adipocytes.


OCR = ETC activity + non-mitochondrial respiration
Best fit between H2 (─) and BMP4 data for direct activation of UCP1

Supplemental Experimental Procedures Basic principles of mathematical modeling in biology
Mathematical modeling is a formalized way of testing mechanistic hypotheses using experimental data and prior knowledge. A mechanistic hypothesis corresponds to an idea of which mechanisms that are essential to produce the observed behavior in the experimental data. For the translation of the hypothesis to such specific models, we use ordinary differential equations (ODEs) of the following format: where represents the states, that corresponds to concentrations of substances, the parameters, that corresponds to the kinetic rate constants, and contains the measured signals, that corresponds to the experimental datasets. The non-linear functions and describe a set of specific dynamic/mechanistic assumptions.
The parameters of biological models have in principle always unknown values (Nyman et al., 2012). We therefore need to define realistic lower and upper limits of their values. Using these limits, we search through a defined space of possible parameters in the search for acceptable parameters. Such acceptable parameters are parameters that give a good enough agreement between model simulations and experimental data. The search for acceptable parameters is an optimization problem based on the least square error: where is the tested values of the parameters, ( ) is referred to as the cost function, ( ) is the measurement data, �( , ) is the simulated curve, ( ) is the standard error of the mean in the measurement data. The summation is over all measured data points. The agreement between model simulation and experimental data can be studied visually and/or formally tested e.g. with a 2 -test (Cedersund and Roll, 2009). We use a 2 -test with 95% confidence to define the set of acceptable parameters.
The outcome of the optimization process decides the next step in the hypothesis testing cycle. If the agreement between model and data is unacceptable from a statistical point of view, the model and corresponding hypothesis are rejected. The rejection is a final conclusion and therefore the next step will be to test another hypothesis.
The other possible outcome from the optimization process is that the agreement between model and data is statistically acceptable, i.e. the model cannot be rejected. This is not a final conclusion, since new datasets may not be in agreement with this so far acceptable model/hypothesis. The next step is therefore to gather all acceptable parameters, i.e. the parameters that give model simulations with a statistical agreement with data, and search for unique predictions that are shared among these parameters. These unique predictions can be used in the design of new experiments, e.g. to be able to discriminate between several acceptable models/hypotheses.

A mathematical model for oxygen consumption in beige and white adipocytes
The mathematical model for oxygen consumption in adipocytes was developed to contain only necessary components to be able to analyze the available data from beige and white adipocytes (Figures 3 and 4).

Basal model equations
The basis of the model is the flow of protons over the inner mitochondrial membrane ( Figure S1). The number of protons in the intermembrane space ( ) is changing with the following ODE: The flow of protons is divided in four parts where the electron transport chain ( ) pumps electrons to increase the proton gradient, and ATPase uses this proton gradient to produce ATP. The two other flows are unspecific leakage ( ) and UCP1-dependent leakage ( 1), also referred to as uncoupled respiration.
The pumping of protons through ETC depends on substrate availability, for example fatty acids (FA), and this effect is potentiated by browning agents (BA1), such as rosiglitazone (Rosi) and BMP4. Several effects of browning agents are allowed in the model (BA1-BA3), and these are handled as continuous parameters optimized for best fit with data. The effects are allowed to be different for different browning agents ( 1 , 2 etc). The flow of protons also depends on the existing proton gradient, which is calculated as the difference between total protons ( ) and . The total proton number depends on added browning agents (BA2), since the total number of mitochondria increases with browning.
The flow of protons down the gradient, used to produce ATP, is assumed to be constant.

=
The leak through the membrane is also assumed to be constant.

=
The uncoupled respiration occurs only if browning agents are added (BA3), and if UCP1 is activated by FA.
BA3 is a continuous parameter optimized for best fit with data and is allowed to be different for different browning agents.
The measured oxygen consumption rate (OCR) is the sum of ETC and non-mitochondrial respiration ( ), i.e. the measured respiration when ETC = 0. The non-mitochondrial respiration depends on the addition of browning agents (BA4).

= ( , ) + • (1 + 4)
Simulating isoproterenol and all-trans retinoic acid stimulation Two different agents are used to stimulate uncoupling: isoproterenol (Iso) that works indirectly via FA, and alltrans retinoic acid (ATRA) that directly binds to UCP1. Both these factors are included in the model as discrete variables that are switched from 0 to 1 when added.
The Iso effect ( ) on FA is described by a single ODE that merges several steps of intracellular signaling from beta-adrenergic receptors to fatty acid release. Two parameters ( and 1 ) were needed for a good agreement between model simulations and data.
The direct effect of ATRA on UCP1 is implemented to be present only if UCP1 is present.
To obtain the observed delayed dynamics for the ATRA response, an ODE for ATRA is used, where the input ( ) is changed from 0 to 1 when ATRA is added and 1 and 2 are optimized for best agreement with data. Two parameters were needed for good agreement with data.

Simulating a mitochondrial stress test
A mitochondrial stress test measures the uncoupled respiration, the maximal respiratory capacity, and the nonmitochondrial respiration by adding agents that block ATPase, increase leakage, and block ETC. In a simulation, these agents are added to the model at the time-points given by the experimental protocol.
In the simulation of a stress test, first oligomycin (OLIGO) is applied to inhibit ATPase.
Because of the dynamics in data, OLIGO is not a simple parameter that is set to 1 when oligomycin is added.
Instead an ODE describes the dynamics of the OLIGO inhibition, where the discrete variable is changed from 0 to 1 when oligomycin is added and the parameter is optimized for best agreement with data. One parameter was enough for good agreement with data.
In the next part of the stress test, DNP is added to increase the unspecific leakage. This value of the input parameter DNP is changed between different simulations to be able to fit with the DNP data (see below section: Model input).

= +
In the final part of the stress test, rotenone (ROT) is added to completely shut down the ETC. The value of ROT is changed from 0 to 1 when added.  The model is simulated to a steady state for the given parameter values and the states are updates to new initial values before each change in the input.

Simulating UCP1 depletion
To be able to fit with the scramble experiments, i.e. the control experiments to the UCP1 depletion measurements in Rosi-treated cells ( Figure S3A-B), an extra parameter was needed ( ). The extra parameter reduces the effect of Iso according to data in Figure 4B in (Bartesaghi et al., 2015). The reason for the reduced effect of Iso can be due to the transfection process.
This extra parameter was estimated for best agreement with data and used only to simulate the UCP1 depletion and scramble experiments ( Figure S3A-B). To simulate UCP1 depletion, 3 is set to 0.01.
The full set of model equations The example values for the BA-parameters show that Rosi has a higher FA-effect at ETC (BA1) than BMP4 and that and BMP4 has a higher increase in proton number and non-mitochondrial respiration (BA2) than Rosi.
The BA3 effect has to do with UCP1 and different values are used for the different tested hypotheses for UCP1 (see section below: Hypothesis testing for BMP4 data).

Model input parameters
Model input parameters are changed between off/on when they are added according to the experimental protocol. The input parameter atra is optimized to fit data. The input parameter DNP will receive different values for different experimental conditions. In white adipocytes (control cells), DNP = 10. In Rosi-treated beige adipocytes, DNP = 5 when iso is present and DNP = 20 when iso is not present. As seen in Figure 3D, the effect of DNP is lower in Rosi-treated cells in the presence of iso. The flexibility of the DNP input parameter will not affect any of the conclusions herein since the conclusions are based on data without DNP (i.e. Figure 4C-E).
The same values for the input parameters are used for both Rosi and BMP4 simulations.
Hypothesis testing for BMP4 data Three hypotheses were tested to explain the BMP4 data: Hypothesis 1 (H1), Hypothesis 2 (H2), and Hypothesis 3 (H3) ( Figure 4B). In the hypothesis testing, parameter values were estimated to be in agreement with control, Rosi, and BMP4 data simultaneously. Also, a qualitative agreement with UCP1 depletion data from Rosi-treated cells should be obtained.