Statistical Physics of Medical Diagnostics: Study of a Probabilistic Model
Alireza Mashaghi Abolfazl Ramezanpour






^aLeiden Academic Centre for Drug Research, Faculty of Mathematics and Natural Sciences, Leiden University, Leiden, The Netherlands
^bHarvard Medical School, Harvard University, Boston, Massachusetts, USA
^cDepartment of Physics, University of Neyshabur, Neyshabur, Iran
*a.mashaghi.tabari@lacdr.leidenuniv.nl




We study a diagnostic strategy which is based on the anticipation of the diagnostic process by simulation of the dynamical process starting from the initial findings. We show that such a strategy could result in more accurate diagnoses compared to a strategy that is solely based on the direct implications of the initial observations. We demonstrate this by employing the mean-field approximation of statistical physics to compute the posterior disease probabilities for a given subset of observed signs (symptoms) in a probabilistic model of signs and diseases. A Monte Carlo optimization algorithm is then used to maximize an objective function of the sequence of observations, which favors the more decisive observations resulting in more polarized disease probabilities. We see how the observed signs change the nature of the macroscopic (Gibbs) states of the sign and disease probability distributions. The structure of these macroscopic states in the configuration space of the variables affects the quality of any approximate inference algorithm (so the diagnostic performance) which tries to estimate the sign/disease marginal probabilities. In particular, we find that the simulation (or extrapolation) of the diagnostic process is helpful when the disease landscape is not trivial and the system undergoes a phase transition to an ordered phase.           



Statistical Physics of Medical Diagnostics: Study of a Probabilistic Model
    Alireza Mashaghi^a,b,*, Abolfazl Ramezanpour^a,c
    December 30, 2023
==========================================================================





§ INTRODUCTION

Statistical physics has been widely used to extract macroscopic properties of a wide range of systems from their microscopic interaction models, yet it has not been employed to medical diagnostics. Given an initial subset of observed signs (symptoms, clinical and laboratory findings) with some prior knowledge about the patient (or a complex system like a biological cell), a diagnosis problem simply asks for the most probable diseases (or macrostates with specific phenotypes) <cit.>. An efficient and accurate diagnostic procedure is important specially in the early stages of diseases, where the number and quality of medical evidences are often insufficient to reach a definite diagnosis. Here, we use approximate inference and optimization algorithms of statistical physics <cit.> to show that a simulation (extrapolation) of the diagnostic process (without doing any real observation) could be helpful as a heuristic strategy in the study of diagnosis problems.


A diagnosis problem usually starts with a (probabilistic) model of (well-defined) sign and disease variables which describes the (statistical) dependencies of the variables; such an effective model of the signs and diseases may come from a microscopic model of the system (human body or biological cell) with emergent macroscopic behaviors that are interpreted as diseases. Various modeling frameworks have been developed and used in medical diagnostics: (i) probabilistic models and belief networks, (ii) neural networks and machine learning methods, and (iii) a complex network approach to the problem.


Bayesian belief networks provide a probabilistic framework to study the sign-disease dependencies <cit.>. The belief networks are represented by tables of conditional probabilities that show how the state of a variable in an acyclic directed graph depends on the state of the parent variables. The above information along with a few simplifying assumptions then are used to infer the marginal sign and disease probabilities for a given set of findings <cit.>. Another approach is to use an artificial neural network to represent the complex relationships of the sign/disease variables <cit.>. The model parameters here are obtained in a learning process using the machine learning techniques <cit.>. Finally, in a network approach to the problem, one constructs a (weighted) symptom-disease network with connections relating the signs to the diseases. This network along with other complementary data, e.g., gene-disease, RNA-disease, protein-disease, metabolite-disease and disease-disease networks, are then utilized as information resources by a diagnostic algorithm <cit.>. 
              
              
In addition to the model, an efficient inference algorithm is needed to estimate the marginal sign and disease probabilities <cit.>. When the number of initially observed signs is too small to make a diagnosis, we need to suggest a number of new medical tests to know the value of several other relevant signs. For this aim, we need an appropriate objective function and optimization algorithm to choose the more informative signs, which can lead us to the right diagnosis with a smaller number of observations. In Refs. <cit.> we proposed probabilistic models of signs and diseases which can systematically take into account the effects of different types of sign-disease, disease-disease, and sign-sign interactions; the models are indeed graphical models of the sign and disease variables with a number of interaction factors, each one connecting a small subset of disease variables to the associated sign variables <cit.>. We introduced approximate inference and optimization algorithms to deal with such probabilistic models, and studied the effects of the model structure and the objective function on the performances of the diagnostic algorithms.

The models we consider are natural generalizations of the simpler probabilistic models studied in previous works <cit.>, which usually assume that only one disease is behind the findings (exclusive diseases assumption), or, the diseases act independently on the signs (causal independence assumption). Moreover, for computational simplicity, it is usually assumed that there is no disease-disease and sign-sign interactions. In Ref. <cit.>, we showed that such interactions can significantly improve the accuracy of diagnosis without resorting to the exclusive diseases or the causal independence assumption.


In this paper, we elaborate more on the nature and behavior of the macroscopic states of the probabilistic models we introduced in the above studies. We employ the mean-field approximation to study the possible changes in the (macroscopic) state of the system as the number of observed signs increases, and to estimate the sign and disease marginal probabilities <cit.>. For the objective function, we choose a function that favors the more polarizing tests, which result in disease probabilities that are closer to zero or one <cit.>. This could be useful especially when the gap between the most probable diseases and the other diseases is small. Moreover, this objective function is easier to compute than a maximum-likelihood function that is typically taken in these problems. Starting from an initial set of observed signs, we use an approximate (Monte Carlo) optimization algorithm to find a sequence of candidate observations that maximizes the above objective function <cit.>. However, instead of the true value of the observed sign at each step, we assume the outcome is given by the most probable value of the sign obtained from the model by the approximate inference algorithm. We show that this strategy is able to improve the quality of diagnosis compared to the case that is merely based on the direct implications of the initial findings. Interestingly, the improvement is observed for nontrivial cases when the system undergoes a phase transition to an ordered phase; i.e., where the effect of observed signs can propagate in the system to influence the state of the other signs (for a similar phenomenon see <cit.>). 





§ MAIN DEFINITIONS

The microscopic state of the system (patient or cell) is identified by the sign values 𝐒={S_i: i=1,…,N_S}, where for simplicity we work with binary sings S_i=± 1. The probability of being in state 𝐒 is given by P(𝐒). The probability P(𝐒) (the model) is parametrized by a set of couplings 𝐊(t), which in general depend on real time t. The conditional probability of the unobserved signs depends on the subset of the observed signs 𝐎(t)={j=1…,N_O(t)} with the values denoted by 𝐒^o(t)={S_j:j ∈𝐎(t)}. 

The macroscopic states (or phenotypes, or diseases) of the system (for large number of signs N_S→∞) can be identified by the Gibbs states of P(𝐒) <cit.>. We label these macroscopic states with 𝐃, for diseases, with 𝐃=0 representing the healthy state. The average of an unobserved sign in state 𝐃 is denoted by ⟨ S_i⟩_𝐃. A pure Gibbs state is characterized by the clustering property; i.e., ⟨ S_iS_j⟩_𝐃-⟨ S_i⟩_𝐃⟨ S_j⟩_𝐃 goes exponentially to zero by the distance of sign nodes i and j in the interaction graph of the sign variables induced by P(𝐒). A mixed state is composed of more than one pure states. In this way, the state of a disease pattern 𝐃 or a cluster of similar disease patterns are represented by the statistical properties of the sign variables in the associated pure or mixed Gibbs states.

We start by asking several interesting questions:



  * What is P(𝐒) and how does it (or the couplings 𝐊) change with time? Here we need a dynamical model of the system to study the stochastic evolution of the sign variables. In the following, we shall assume some reasonable structures for the model and leave this problem for future studies. Instead of going from the model P(𝐒) to the macroscopic states, we start from the diseases 𝐃 and obtain the model from the joint probability of the sign and disease variables P(𝐒;𝐃)=P(𝐒|𝐃)P_0(𝐃). Then, the model is obtained by summing over the disease variables P(𝐒)=∑_𝐃P(𝐒|𝐃)P_0(𝐃). The conditional probability P(𝐒|𝐃) can be a decreasing function of the distance of 𝐒 and a reference sign configuration 𝐒(𝐃). Here the S_i(𝐃) represent the most probable symptoms of disease 𝐃. These models could be useful (in the absence of the realistic models) as benchmarks in the study of a diagnosis problem. 
   

  * Do we see a significant change of behavior with time in P(𝐒)? For example, from weak sign correlations to a regime of strong correlations. Typically, we encounter strong correlations close to a phase transition from one macroscopic state to another state. As we will see, even simple (but plausible) models of signs and diseases can exhibit both continuous and discontinuous phase transitions as the strength of the sign and disease interactions are varied. In particular, the phase transition can be induced by increasing the number of observed signs for given strength of the interactions.   


  * How does the structure of P(𝐒) affect the diagnosis? Here we need an efficient and approximate inference algorithm to compute the sign and disease probabilities. It is easy to obtain very good estimations of these marginal probabilities as long as there is only one macroscopic (pure) state, or there are a number of symmetry-related states. Otherwise, the above algorithms will not converge or will need a very large computation time to provide a fair sampling of the probability distribution. We will see how the convergence and quality of an approximate inference algorithm which is based on the mean-field approximation affect the diagnostic performance. 



Consider a set of N_D binary variables 𝐃={D_a=0,1: a=1,…,N_D}, where D_a=0,1 shows the absence or presence of disease a. Each disease is assigned a positive weight W_a, to take into account the significance of diseases. In the following we assume the W_a are uniformly distributed in (0,1). The joint probability distribution of the sign and disease variables (i.e., the model) is identified by P(𝐒;𝐃)=P(𝐒|𝐃)P_0(𝐃). Here P_0(𝐃) is the prior probability distribution of diseases, which could depend on the patient's characteristics such as gender and age and disease properties such as duration of a disease, mortality rate and transmission rate among others. In the following, we shall assume the prior probability is factorized as P_0(𝐃) ∝exp(∑_a K_a^0D_a). The parameters K_a^0 can also be used to control the expected number of present diseases.

Let P_true(𝐒|𝐃) to be the true (or empirically estimated) probability distribution of the sign variables given the disease hypothesis 𝐃. In practice, we may have access only to a small subset of marginal probabilities of the above probability distribution. For instance, suppose we are given the sign probabilities P_true(S_i|nodisease), P_true(S_i,S_j|onlyD_a), and P_true(S_i,S_j|onlyD_a,D_b) conditioned on the presence of no diseases, the presence of only one disease, and the presence of only two diseases, respectively. Using the maximum entropy principle, for the conditional probability distribution of the signs we take

    P(𝐒|𝐃)=1/Z(𝐃)ϕ_0(𝐒) ×∏_a ϕ_a(𝐒|D_a)×∏_a < bϕ_ab(𝐒|D_a,D_b),

where the partition function Z(𝐃) is obtained from normalization. 

The disease interaction factors (ϕ_0, ϕ_a, ϕ_ab) can in general be parametrized by the couplings of all the possible multi-sign interactions. As customary in maximum entropy modeling, assuming an exponential family, the parameters
sufficient to describe the above marginal probabilities are involved in the one-sign terms (K_i^0S_i), the one-disease-one-sign interactions (K_i^aD_aS_i), the one-disease-two-sign interactions (K_ij^aD_aS_iS_j), the two-disease-one-sign interactions (K_i^abD_aD_bS_i), and finally the two-disease-two-sign interactions (K_ij^abD_aD_bS_iS_j). More precisely, the disease interaction factors are given by 

    ϕ_0(𝐒)    = e^∑_i K_i^0 S_i,
    ϕ_a(𝐒|D_a)    = e^D_a[∑_i K_i^a S_i+∑_i<j K_ij^a S_iS_j],
    ϕ_ab(𝐒|D_a,D_b)    = e^D_aD_b[∑_i K_i^ab S_i+∑_i<j K_ij^ab S_iS_j].

Figure <ref> shows the interaction graph of the sign and disease variables related by the above interaction factors. We use M_a,k_a and M_ab,k_ab for the number and connectivity of one-disease and two-disease interaction factors, respectively.

In principle, the information provided by the marginal probabilities of the true (or empirical) probability distribution is sufficient to determine the model parameters <cit.>. 
Note that ϕ_0(𝐒) is responsible for the probability of observing 𝐒 in the absence of any diseases, where the most probable value is S_i=-1. It is reasonable to assume that in the healthy case each sign takes the positive value with a small probability independently of the other sign values.  

The joint probability distribution of the sign and disease variables can be rewritten as, P(𝐒;𝐃) ∝exp(-ℋ(𝐒;𝐃)), where the energy function reads as follows

    ℋ(𝐒;𝐃) = -∑_i h_i(𝐃) S_i-∑_i<jJ_ij(𝐃)S_iS_j
    +ln Z(𝐃)-ln P_0(𝐃).

Here, the partition function and the new couplings are:

    Z(𝐃)    =∑_𝐒 e^∑_i h_i(𝐃) S_i+∑_i<jJ_ij(𝐃)S_iS_j,
    
    h_i(𝐃)    = K_i^0+∑_aK_i^aD_a+∑_a<bK_i^abD_aD_b,
    
    J_ij(𝐃)    = ∑_aK_ij^aD_a+∑_a<bK_ij^abD_aD_b.

From the above model, we can extract simpler models depending one the maximum number of disease and sign variables that are involved in the interactions; for instance, we could have the D1S1 (one-disease-one-sign), D1S2 (one-disease-two-sign), D2S1 (two-disease-one-sign), and D2S2 (two-disease-two-sign) models.

In the following, we consider only the D1S1 and D2S1 models, where we can exactly compute the partition function 
Z(𝐃)=∏_i ( 2cosh h_i(𝐃)).
For these models, we can also exactly compute the model parameters given the true marginal probabilities,

    K_i^0    =1/2ln(P_true(S_i=+1|nodisease)/P_true(S_i=-1|nodisease)),
    
    K_i^a   =1/2ln(P_true(S_i=+1|only D_a)/P_true(S_i=-1|only D_a))-K_i^0,
    
    K_i^ab   =1/2ln(P_true(S_i=+1|only D_a,D_b)/P_true(S_i=-1|only D_a,D_b))-K_i^0-K_i^a-K_i^b.
  


For a given subset 𝐎 of observed signs with values 𝐒^o, the disease probabilities are obtained from 

    P(D_a=1|𝐒^o)=1/𝒵(𝐒^o)∑_𝐃 D_a e^-ℋ(𝐃|𝐒^o),

where ℋ(𝐃|𝐒^o) = -logℒ(𝐃|𝐒^o) is the log-likelihood function

    ℋ(𝐃|𝐒^o)=  -∑_a K_a^0D_a-
    ∑_i∈𝐎S_i^oh_i(𝐃)+∑_i∈𝐎ln(2cosh h_i(𝐃) ),

and 𝒵(𝐒^o)= ∑_𝐃exp(-ℋ(𝐃|𝐒^o)) is another normalization constant.
As before, the prior probability distribution is P_0(𝐃) ∝exp(∑_a K_a^0 D_a).
It is easy to show that the marginal probability of an unobserved sign is given by:

    P(S_i=1|𝐒^o)=1/𝒵(𝐒^o)∑_𝐃(1+tanh h_i(𝐃)/2) e^-ℋ(𝐃|𝐒^o).

The approximate equations for the D1S2 and D2S2 models can be found in <cit.>.


 






§ THE HOMOGENEOUS FULLY-CONNECTED MODELS

The direct problem of inferring the marginal sign/disease probabilities from the above models can be solved exactly as long as the model parameters do not depend on the sign or disease labels.
The thermodynamic limit here is defined by the limit N_D,N_S,N_O →∞ such that γ= N_D/N_S and n_o=N_O/N_S remain finite. To provide some order of magnitude, it is useful to mention that in Internist (a probabilistic model for internal diseases <cit.>) the number of diseases is about 500 and the number of associated signs is around 4000. In addition, we need to scale the model parameters as K_ij^a=κ_ij^a/(N_SN_D), K_ij^ab=κ_ij^ab/(N_SN_D^2), K_i^0=κ_i^0, K_i^a=κ_i^a/N_D, K_i^ab=κ_i^ab/N_D^2, and K_a^0=κ_a^0; the scaling ensures that the energy function is extensive (proportional to N_S).


The sign and disease probabilities are obtained by minimizing the following free energy with respect to x=P(D=1) and y_u=P(S=1) (for an unobserved sign), 

    f(x,y)=-γ𝒮(x)-(1-n_o)𝒮(1+y_u/2)+𝒮(1+z(x)/2)
    
    -h(x)(y-z(x))-1/2J(x)(y^2-z^2(x))-γκ_a^0 x.

The value of the observed signs enters in y=n_o y_o+(1-n_o)y_u with y_o=(∑_i∈𝐎S_i^o)/N_O, and y_u=(∑_i∉𝐎S_i)/(N_S-N_O). Here z is the solution to z=tanh(h(x)+J(x) z), and 𝒮(x)=-xlog x-(1-x)log(1-x) is the Gibbs-Shannon entropy function. Moreover, the effective field h(x) =κ_i^0+κ_i^a x+1/2κ_i^abx^2 and the coupling J(x) =κ_ij^a x+1/2κ_ij^abx^2 (see Appendix <ref> for the derivations). Each local or global minimum of the free energy can be considered as a macroscopic state of the system.  Figure <ref> shows how the sign and disease probabilities change with the number of observations, when all the observed signs have a positive value (see also Fig. <ref> in Appendix <ref>). As the figures show, a new macroscopic state can appear continuously or discontinuously depending on the value of the model parameters.




§ THE INHOMOGENEOUS MODELS: MEAN-FIELD APPROXIMATION

In this section, we find an estimation of the sign and disease probabilities for arbitrary couplings 𝐊. To this end, we write D_a=⟨ D_a ⟩+δ D_a and S_i=⟨ S_i ⟩+δ S_i where the δ D_a=D_a-⟨ D_a ⟩ and δ S_i=S_i-⟨ S_i ⟩ are small deviations from the mean values. The mean-field (MF) approximation here is obtained by neglecting the second order deviations in a Taylor expansion around the mean values <cit.>. In the following, we shall restrict ourselves to the D1S1 and D2S1 models, where the normalization function Z(𝐃) can be computed exactly; for the D1S2 and D2S2 models we need also to compute this function within the MF approximation (see Appendix <ref>).

In this way, the MF approximation for the sign and disease probabilities are obtained by solving the self-consistency equations x_a =exp(h_a(𝐱))/(1+exp(h_a(𝐱))), with P(D_a=1)=x_a and P(S_i=1)=(1+tanh h_i(𝐱))/2. Here, the effective fields experienced by the sign and disease variables are given by

    h_i(𝐱)    = K_i^0+ ∑_a K_i^a x_a+∑_a<bK_i^abx_ax_b,
    
    h_a(𝐱)    = K_a^0+∑_i∈𝐎[S_i^o-tanh(h_i(𝐱))](K_i^a+∑_b≠ aK_i^abx_b).

The equations are solved by iteration starting from random initial values for the x_a. The time complexity of this algorithm is of order N_ON_D^3 in a fully-connected model. The fixed points of these equations are considered as the macroscopic states of the system (Gibbs states). As long as there is only one macroscopic state, the iteration algorithm converges easily to the single fixed point of the equations. Non-convergence of the iteration algorithm is a signature of the presence of more than one fixed point.

To check the performances of the algorithms, we shall assume that the true model is given by an exponential probability distribution P_true(𝐒|𝐃) ∝exp(-β H(𝐒;𝐒(𝐃))). Here 𝐒(𝐃) gives the most probable symptoms of disease pattern 𝐃, and H(𝐒;𝐒(𝐃)) is the Hamming distance (number of different elements) of the two sign configurations. Moreover, β is a positive parameter that controls the structure of the true model around the symptoms 𝐒(𝐃); the diseases are more clearly distinguished for larger values of β. We assume that each element S_i(𝐃) (for i=1…,N_S) takes the positive and negative values with equal probability, except for the healthy case (𝐃=0) where all the elements are negative. Given the true model, we use the true marginal probabilities to construct e.g. the D2S1 model. 

Suppose that we are given a subset 𝐎 of N_O observed signs with values 𝐒^o. A simple diagnostic procedure works by computing the posterior disease probabilities conditioned on the observations P(D_a|𝐒^o). Then, the most probable diseases or those that have a probability greater than a threshold value, are reported as the diagnosed diseases; in the following, we shall assume that the most probable diseases, within a small window of size δ P_D=0.01, are the present ones. Figure <ref> displays the accuracy of such a diagnosis with the D1S1 and D2S1 models for a small number of sign and disease variables. The figure also shows the probability gap between the most probable disease(s) and the other diseases. A patient with disease pattern 𝐃 and N_O initial observed signs from the most probable symptoms 𝐒(𝐃) is presented to the model for diagnosis; a disease pattern is chosen with a probability proportional to the weights W_a of the present diseases in 𝐃. From <cit.> we know that the D1S1 and D2S1 models work well so long as the number of present diseases in 𝐃, denoted by |𝐃|, is less than or equal to two; that is why we choose patients with a small number of diseases. As the figure shows, we obtain more accurate predictions as the parameter β increases. The situation is different when we have to resort to an approximate inference algorithm. We see in Fig. <ref> that the MF approximation does not provide accurate estimations of the sign and disease marginal probabilities for large β, where the algorithm does not converge. Here the best performances are observed for intermediate values of β. 

    


§ DIAGNOSIS BY SIMULATION OF THE DIAGNOSTIC PROCESS

It may happen that the information provided by the initial number of observations are not enough to reach a reliable diagnosis, especially in the early stages of the diseases. Thus, we need a good strategy to choose the most informative signs for the next observations. Here the goal could be to reach the right diagnosis with a minimum number of the medical tests <cit.>. Thus, for the objective function we propose an increasing function of the polarizations (deviations from the neutral value) in the posterior disease probabilities. The optimal choice then is provided by the most polarizing observation conditioned on the value of the previous observations. In contrary to the maximum likelihood function which is computationally hard to compute, the above objective function can easily be computed given the posterior disease probabilities. And, one can easily incorporate the importances of the diseases (the W_a) into the objective function, to assign more weight to polarization of the more important diseases. 

In a sequential diagnostic process of length T, we do the medical tests one by one and at each step we obtain the true value of the observed sign (this is called Diags-I in <cit.>). To obtain an optimal sequence of medical tests, one has to simulate in advance a diagnostic process of T observations without doing any real observation (this is called Diags-II in <cit.>). This simulation of the observation process, or extrapolation from the initial set of observations, is proposed here as another heuristic approach to diagnosis to fully exploit the statistical dependencies of the variables provided by the model in addition to the initial medical tests. To this end, we use the mean-field approximation to compute the posterior disease probabilities. Then, the (Monte Carlo) optimization algorithm of Ref. <cit.> is used to maximize an objective functional of T observations, 

    ℰ[𝐎(T)]= ∑_t=1^T (1/∑_a W_a∑_a W_a | P(D_a=1)-1/2| ).

Here 𝐎(T)={j_1,…,j_T} is the sequence of observations. One can also add the cost or relevance of the observed signs to this objective function <cit.>. The optimization algorithm starts from a random sequence of T observations, uses the marginal sign probabilities to generate a sequence of new observations, and accepts the suggested sequence if the objective function increases.    

Note that the above problem is indeed a stochastic optimization problem, where the objective function depends on the stochastic outcomes of the observations <cit.>. To simplify the computation, we assume that the observed sign j at each time step takes the most probable value identified by the marginal probability P(S_j|𝐒^o) conditioned on the value of the previous observed signs.  

Figure <ref> shows how the above objective function and the optimization algorithm perform. The figure displays the changes the first right and wrong diagnosis times compared to a random sequence of T observations  <cit.>; the first right diagnosis time T_R is the first time (number of observations) the probability of having a right disease becomes larger than a threshold value, here P_th=0.9. Similarly we define the first wrong diagnosis time T_W. In Fig. <ref>, we compare the accuracy of the diagnosis with the D2S1 model before and after extrapolation for T=N_O/2 steps. Here the sign/disease marginal probabilities are computed by the MF approximation. Similar comparisons are shown also in Figs. <ref> and <ref>.






§ DISCUSSION

In summary, depending on the model and the strength of the model parameters, new macroscopic states can appear as the number of observed signs increases. This could be helpful because the disease probabilities are usually more informative within such states. On the other hand, this affects the algorithm convergence and consequently the quality of the sign and disease probabilities which are computed by the approximate inference algorithm. More advanced and accurate algorithms can of course improve the quality of inference, but at the expense of more computational time <cit.>.    

We showed that simulation of the diagnostic process provides a useful strategy for diagnosis when a naive approach that is based on the direct implications of the observed signs is not very helpful. In other words, this strategy works in the ordered phase of the system where the values of the observed signs significantly affect the probability distribution of the unobserved signs; the classical example is a ferromagnetic spin system in the ordered (low temperature) phase where the values of the boundary spins determine the physical (Gibbs) state of the system. In this way, we can define a critical number of initial observations which are needed to enter such an ordered state, for systems that display a phase transition.   

Here, for the sake of efficiency, we assumed that each "observed" sign in the simulation takes the most probable value predicted by the model. Moreover, we used a very naive optimization algorithm to find the optimal sequence of observations. A more accurate study should consider the stochastic nature of the simulated observations, and employ a more sophisticated optimization algorithm, e.g., simulated annealing. Finally, it would be interesting to have a microscopic (or phenomenological) model of patient (or an ensemble of patients) to study the time evolution of the sign probability distribution, and the emergent macroscopic (disease) states.           



§ PERSPECTIVES

An accurate medical diagnosis from a limited number of findings (e.g. at the early stages of diseases) should exploit all the statistical information on the sign/disease dependencies observed in the clinical and laboratory data. Such interdependencies are emerging due to the advancements in omics technologies and progress in population studies and aging research (e.g. identification of co-occurrence of age-related diseases). We note that the existing datasets lack the necessary probabilistic information needed for our approach, as such new data need to be generated. These studies will be the subject of our future works and in the current article, we are primarily addressing the mathematics and statistical physics communities.

Let us recall briefly the kind of statistical data we need to construct the models studied in this paper. First, note that these models have been obtained from an expansion around the healthy state where the number of involved diseases is small (|𝐃|=1,2) <cit.>. On the other hand, given a disease hypothesis 𝐃, it is usually assumed that the sign variables are uncorrelated in a zero-order approximation of the signs statistics <cit.>. Here, the necessary data are encoded in the conditional probabilities P(S_i|onlyD_a) (in D1S1 model) or P(S_i|onlyD_a,D_b) (in D2S1 model). Obviously, we expect to have two-sign correlations, or higher-order sign correlations, even in presence of only a single disease. But taking into account these correlations considerably increases the computational complexity of the problem. Additionally, it is in practice very difficult to obtain statistically good clinical data which capture the higher-order correlations. Nevertheless, in the end, it is the collection of available empirical data that determines the structure of the model. 

The method can in principle be applied to any diagnostic problem to infer the macroscopic state (phenotype) of the system from a limited number of evidences. This could be, for instance, the problem of assigning a state to a biological cell or a complex electronic device. In particular, assignment of state to a cell is a major challenge in immunology and cancer biology and it has complicated developing therapies for cancer and autoimmunity. We envision that our approach will be generically applied to a wide range of problems in medicine, science and technology.













§ THE HOMOGENEOUS FULLY-CONNECTED MODELS

As long as the model parameters do not change with the sign or disease labels, we can write all the quantities in terms of the collective variables x=(∑_aD_a)/N_D and y=(∑_iS_i)/N_S. Then for large number of signs (N_S→∞), we get

    1/N_Sln Z(𝐃)    ≈𝒮(1+z(x)/2)+h(x)z(x)+1/2J(x)z^2(x),
    1/N_Sℋ(𝐒;𝐃)    ≈ -h(x)(y-z(x))-1/2J(x)(y^2-z^2(x))+𝒮(1+z(x)/2)-γκ_a^0 x,

where γ=N_D/N_S, and

    h(x)    =κ_i^0+κ_i^a x+1/2κ_i^abx^2,
    
    J(x)    =κ_ij^a x+1/2κ_ij^abx^2.

Here we take the scaling 

    K_i^0    =κ_i^0,0.5cm  K_a^0=κ_a^0,
    
    K_i^a    =1/N_Dκ_i^a,0.5cm K_i^ab=1/N_D^2κ_i^ab,
    
    K_ij^a    =1/N_SN_Dκ_ij^a,0.5cm K_ij^ab=1/N_SN_D^2κ_ij^ab.

Moreover, z is the solution to

    z=tanh(h(x)+J(x) z).

which minimizes the following free energy

    f(z) = -𝒮(1+z/2)-h(x)z-1/2J(x)z^2,
  
Here, for brevity, we defined the Shanon entropy function

    𝒮(p)=-pln p-(1-p)ln(1-p).



To take into account the value of the observed signs, we write y=n_o y_o+(1-n_o)y_u with y_o=(∑_i∈𝐎S_i^o)/N_O, y_u=(∑_i∉𝐎S_i)/(N_S-N_O), and n_o=N_O/N_S. In this way, the grand partition function is given by

    𝒵(𝐒^o) ≃∫_0^1 dx ∫_-1^1 dy_u e^-N_S f(x,y),

At the end, the self-consistency equations for the x and y_u variables in the thermodynamic limit (N_S→∞), are obtained by minimizing the following free energy 

    f(x,y)=-γ𝒮(x)-(1-n_o)𝒮(1+y_u/2)
    
    -h(x)(y-z(x))-1/2J(x)(y^2-z^2(x))+𝒮(1+z(x)/2)-γκ_a^0 x.

Figure <ref> shows how the above free energy behaves when the model parameters in the D2S1 model are varied.




§ THE INHOMEGEUOUS MODELS: MEAN-FIELD APPROXIMATION

For the D1S1 and D2S1 models we can compute some quantities exactly, therefore, we present the mean-field approximation for these models separately.
 


 §.§ In the absence of the sign-sign interactions

Here we write D_a=x_a+δ D_a where x_a=⟨ D_a ⟩, and δ D_a=D_a-x_a is a small deviation from the mean value.
Then, the local field experienced by sign i is

    h_i(𝐃)= h_i(𝐱)+ ∑_a K_i^a δ D_a+∑_a<bK_i^ab(δ D_ax_b+x_aδ D_b+δ D_aδ D_b)=h_i(𝐱)+δ h_i,

where

    h_i(𝐱)    = K_i^0+ ∑_a K_i^a x_a+∑_a<bK_i^abx_ax_b.

The Hamiltonian can be written as

    ℋ(𝐃|𝐒^o) =  -∑_a K_a^0D_a-
    ∑_i∈𝐎S_i^o[h_i(𝐱)+δ h_i]+∑_i∈𝐎ln(2cosh [h_i(𝐱)+δ h_i] ),

Expanding the last term up to the first order deviations δ D_a, we get

    ℋ(𝐃|𝐒^o) ≈ℋ_0-∑_a h_a(𝐱) D_a,

with

    h_a(𝐱)    = K_a^0+∑_i∈𝐎[S_i^o-tanh(h_i(𝐱))](K_i^a+∑_b≠ aK_i^abx_b).


Then, the average values x_a are obtained by the following self-consistency equations:

    x_a =e^h_a(𝐱)/1+e^h_a(𝐱),

The equations are solved by iteration starting from random initial values for the x_a. 




 §.§ In the presence of the sign-sign interactions

In general, the partition function in the MF approximation reads

    Z(𝐃) ∝∏_i (2cosh[h_i(𝐃)+∑_j≠ iJ_ij(𝐃)z_j(𝐃)] ),

where the z_i are solutions to

    z_i(𝐃)=tanh[h_i(𝐃)+∑_j≠ iJ_ij(𝐃)z_j(𝐃)].

Define δ D_a=D_a-⟨ D_a ⟩ and δ S_i=S_i-⟨ S_i ⟩. For brevity, we take x_a=⟨ D_a ⟩ and y_i=⟨ S_i ⟩.
Note that y_i=S_i^o is fixed for i∈𝐎 (the subset of observed signs).

Then, to first order in the δ D_a and δ S_i, we have

    ℋ(𝐒;𝐃) ≈ℋ_0 -∑_i∉𝐎 [h_i(𝐱)+∑_j≠ iJ_ij(𝐱)y_j] S_i-∑_a h_a(𝐱,𝐲,𝐳) D_a,

where

    h_a(𝐱,𝐲,𝐳)=K_a^0+∑_i B_i^a(𝐱) y_i+∑_i<j B_ij^a(𝐱) y_iy_j 
    
    -∑_itanh[h_i(𝐱)+∑_j≠ iJ_ij(𝐱)z_j]
    ×( B_i^a(𝐱)+∑_j≠ i[B_ij^a(𝐱)z_j+J_ij(𝐱)χ_j^a])

Here, the new introduced local fields are

    B_i^a(𝐱)    =K_i^a+∑_b≠ aK_i^abx_b,
    
    B_ij^a(𝐱)    =K_ij^a+∑_b≠ aK_ij^abx_b,

and the susceptibility χ_i^a is given by 

    χ_i^a = ∂ z_i/∂ x_a=(1+tanh^2[h_i(𝐱)+∑_j≠ iJ_ij(𝐱)z_j])
    ×( B_i^a(𝐱)+∑_j≠ i[B_ij^a(𝐱)z_j+J_ij(𝐱)χ_j^a]).


In summary, the mean-field equations read as follows,

    x_a    =e^h_a(𝐱,𝐲,𝐳)/1+e^h_a(𝐱,𝐲,𝐳),
    
    y_i    =S_i^o,0.5cm i∈𝐎
    
    y_i    =tanh[h_i(𝐱)+∑_j≠ iJ_ij(𝐱)y_j],0.5cm i∉𝐎
    
    z_i    =tanh[h_i(𝐱)+∑_j≠ iJ_ij(𝐱)z_j].

We solve the equations by iteration starting from random initial values for the x_a, y_i (i∈𝐎),z_i, and the χ_i^a. 



prsty


Ledley Ledley, R.S., Lusted, L.B. (1959) Reasoning foundations of medical diagnosis; symbolic logic, probability, and value theory aid our understanding of how physicians reason. Science 130:9-21.

MillerDDSS Miller, R.A. and Geissbuhler, A. (1999). Clinical Diagnostic Decision Support Systems-An Overview, Springer New York, Page 3-34, ISBN: 978-1-4757-3903-9.

Current-2016 Papadakis, M., McPhee, S.J., Rabow, M.W. (2016). Current Medical Diagnosis and Treatment. 55 edition, LANGE CURRENT Series, 1920 pages, ISBN: 0071845097.




HR-book-2006Hartmann, Alexander K., and Heiko Rieger, eds. New optimization algorithms in physics. John Wiley & Sons, 2006.

MM-book-2009 Mezard, M., and Montanari, A. (2009). Information, physics, and computation. Oxford University Press.



Spiegelhalter Spielgelharter, D. J. (1987). Probabilistic Expert Systems in Medicine. Statistical Science, 2, 3-44.

internist-1 Shwe, M. A., Middleton, B., Heckerman, D. E., Henrion, M., Horvitz, E. J., Lehmann, H. P., and Cooper, G. F. (1991). Probabilistic diagnosis using a reformulation of the INTERNIST-1/QMR knowledge base. Methods of information in Medicine, 30(4), 241-255.

Heckerman-1 Heckerman, D. E., and Shortliffe, E. H. (1992). From certainty factors to belief networks. Artificial Intelligence in Medicine, 4(1), 35-52.

Miller94 Miller, R.A. (1994). Medical diagnostic decision support systems—past, present, and
future: a threaded bibliography and commentary. J Am Med Inform Assoc 1:8–27.

Nikovski-ieee-2000 Nikovski, D. (2000). Constructing Bayesian networks for medical diagnosis from incomplete and partially correct statistics. Knowledge and Data Engineering, IEEE Transactions on, 12(4), 509-516.



B-nc-1990 Baxt, William G. "Use of an artificial neural network for data analysis in clinical decision-making: the diagnosis of acute coronary occlusion." Neural computation 2.4 (1990): 480-489.

P-ieee-1998 Penedo, Manuel G., et al. "Computer-aided diagnosis: a neural-network-based approach to lung nodule detection." IEEE Transactions on Medical Imaging 17.6 (1998): 872-880.

Khan-nm-2001 Khan, Javed, et al. "Classification and diagnostic prediction of cancers using gene expression profiling and artificial neural networks." Nature medicine 7.6 (2001): 673.



Murphy-book Murphy, K. P. (2012). Machine learning: a probabilistic perspective. MIT press.



Goh-pnas-2007 Goh, Kwang-Il, et al. "The human disease network." Proceedings of the National Academy of Sciences 104.21 (2007): 8685-8690.

Barabasi-nr-2011 Barabási, Albert-László, Natali Gulbahce, and Joseph Loscalzo. "Network medicine: a network-based approach to human disease." Nature reviews. Genetics 12.1 (2011): 56.


G-gm-2014 Gustafsson, Mika, et al. "Modules, networks and systems medicine for understanding disease and aiding diagnosis." Genome medicine 6.10 (2014): 82.

Sun-bmc-2014 Sun, Kai, et al. "Predicting disease associations via biological network analysis." BMC bioinformatics 15.1 (2014): 304.

Liu-sr-2015 Liu, Wei, et al. "Integrative analysis of human protein, function and disease networks." Scientific reports 5 (2015).

Suratanee-2015 Suratanee, A., and Plaimas, K. (2015). DDA: A Novel Network-Based Scoring Method to Identify Disease–Disease Associations. Bioinform Biol Insights. 9: 175–186.


cooper-ai-1990 Cooper, G. F. (1990). The computational complexity of probabilistic inference using Bayesian belief networks. Artificial intelligence, 42(2), 393-405.

Jordan Jordan, M. I. (2004). Graphical models. Statistical Science, 140-155.



paper-1 Abolfazl Ramezanpour and Alireza Mashaghi. Toward First Principle Medical Diagnostics: On the Importance of Disease-Disease and Sign-Sign Interactions. Frontiers in Physics 5: 32 (2017).

paper-2 Abolfazl Ramezanpour and Alireza Mashaghi. Uncovering a hidden disease pattern by simulating the clinical diagnostic process. Scientific Reports 8(1):2436 (2018).


MF-book-2001 Opper, Manfred, and David Saad, eds. Advanced mean field methods: Theory and practice. MIT press, 2001.

MM-jstat-2006 Mézard, Marc, and Andrea Montanari. "Reconstruction on trees and spin glass transition." Journal of statistical physics 124.6 (2006): 1317-1350.


gibbs-pnas-2007 Krzakała, Florent, et al. "Gibbs states and the set of solutions of random constraint satisfaction problems." Proceedings of the National Academy of Sciences 104.25 (2007): 10318-10323.



KR-nc-1998 Kappen, H. J., and Rodriguez, F. B. (1998). Efficient learning in Boltzmann machines using linear response theory. Neural Computation, 10(5), 1137-1156.

T-pre-1998 Tanaka, T. (1998). Mean-field theory of Boltzmann machine learning. Physical Review E, 58(2), 2302.

T-jstat-2012 Ricci-Tersenghi, F. (2012). The Bethe approximation for solving the inverse Ising problem: a comparison with other inference methods. Journal of Statistical Mechanics: Theory and Experiment, 2012(08), P08015.

inverse-2017 Nguyen, H. C., Zecchina, R., Berg, J. Inverse statistical problems: from the inverse Ising
problem to data science. Advances in Physics, 66 (3), 197-261 (2017).



BL-book-1997 Birge, J. R., and Louveaux, F. (2011). Introduction to stochastic programming. Springer Science & Business Media.

ABRZ-prl-2011 Altarelli, F., Braunstein, A., Ramezanpour, A., and Zecchina, R. (2011). Stochastic matching problem. Physical review letters, 106(19), 190601.


































