Below is a step-by-step outline showing how one can:
	1.	Derive (or at least set up) the apparent solubility vs. pH expression for a tetra-protic molecule (i.e., having 4 labile protons) in the presence of a strong base.
	2.	Incorporate known/unknown parameters (one neutral solubility, three solubility products for 1:1, 1:2, 1:3 salts, and four pKa values).
	3.	Write a Python program that uses these equations to fit experimental solubility vs. pH data, thus determining (or refining) the pKa’s, solubilities, and solubility-product constants.

1. System Description

Let the compound be , having four acidic protons with known (or to-be-fitted) pKa values:

\[
\begin{aligned}
\mathrm{H_4A}
&\xrightleftharpoons{\mathrm{p}K_{a1}} \mathrm{H_3A}^{-} + \mathrm{H}^+ \\
\mathrm{H_3A}^{-}
&\xrightleftharpoons{\mathrm{p}K_{a2}} \mathrm{H_2A}^{2-} + \mathrm{H}^+ \\
\mathrm{H_2A}^{2-}
&\xrightleftharpoons{\mathrm{p}K_{a3}} \mathrm{HA}^{3-} + \mathrm{H}^+ \\
\mathrm{HA}^{3-}
&\xrightleftharpoons{\mathrm{p}K_{a4}} \mathrm{A}^{4-} + \mathrm{H}^+
\end{aligned}
\]

We also assume:
	1.	Neutral solid:  with intrinsic solubility .
	2.	Salt forms precipitate with a strong base cation (, e.g. ). We consider three salt forms whose solubility products are known (or to be fitted):
	•	 with 
	•	 with 
	•	 with 

(If there were a fully deprotonated salt  with a known , you could similarly incorporate it. In this problem, we assume either it is too soluble to precipitate or is out of the pH range of interest.)

2. Equilibria and pH Dependence

2.1. Dissociation Fractions

For a generic tetra-protic acid , in solution at a given pH (with ), the ratio of each species can be written in terms of the four pKa values.  One often defines (in the absence of precipitation) the “alpha” fractions:



Here,
	•	 is the fraction of 
	•	 is the fraction of 
	•	 is the fraction of 
	•	 is the fraction of 
	•	 is the fraction of .

However, in real solubility calculations, once any solid phase precipitates, the solution concentrations of these species become constrained by whichever solid phase is controlling solubility at that pH.

2.2. Phases That Can Precipitate

We have potentially four solids:
	1.	Neutral form , with an intrinsic solubility .
	•	If this alone were controlling solubility, we would have:

and the total dissolved concentration would be
\[
S_{\mathrm{neutral}}(H)
\;=\; S_0 \bigl(\alpha_0 + \alpha_1 + \alpha_2 + \alpha_3 + \alpha_4 \bigr)
\quad \text{(computed at equilibrium)}
\]
with the appropriate expression for the  ratios.
	2.	1:1 salt  with .
	•	If this salt is controlling solubility, then at equilibrium:

Given the pH (thus known ) and the strong base  (some function of total added base, though in many models we approximate \([\mathrm{M^+}}] \approx\) constant or at least known from charge balance), one obtains .  Then the other species  follow from the acid-base equilibria.  Summing these yields

	3.	1:2 salt  with .
	•	At saturation:

From that, again use the pKa equilibria to find  and get

	4.	1:3 salt  with .
	•	Similarly:

Then get the total dissolved -species from all equilibria.

(If needed, one could also include the 1:4 salt, .)

2.3. Overall Apparent Solubility

Because any one of these solid phases (neutral or salt) can “take over” as the limiting (lowest) solubility phase at a given pH, the actual or apparent solubility  at pH  pH is:

\[
S_{\mathrm{app}}(\mathrm{pH})
\;=\;
\min\Bigl[
S_{\mathrm{neutral}}(H),\;
S_{\mathrm{1{:}1\,salt}}(H),\;
S_{\mathrm{1{:}2\,salt}}(H),\;
S_{\mathrm{1{:}3\,salt}}(H)
\Bigr].
\]

This leads to a piecewise curve across pH.  In many practical situations:
	•	At low pH, the neutral form  controls (lowest solubility).
	•	At some intermediate pH, the 1:1 salt takes over.
	•	Then possibly the 1:2 salt at even higher pH, etc.

Hence, plotting  vs. pH often results in multiple “legs” each controlled by a different phase.

3. Fitting Experimental Data

Suppose you have experimental data , a set of pH values and measured total solubilities.  Your unknowns might be:
	•	 (intrinsic neutral solubility)
	•	
	•	

You can perform a non-linear least-squares fit of your theoretical curve



to the observed data .

3.1. Implementation Outline
	1.	Define a Python function app_solubility(pH, params) that:
	•	Extracts parameters .
	•	Computes .
	•	For each possible controlling phase (neutral, 1:1 salt, 1:2 salt, 1:3 salt):
	1.	Compute the saturating concentration of that key species (e.g.,  for neutral, or , etc.).
	2.	Use the acid-base equilibrium relationships to convert from that “key species concentration” to total dissolved .
	•	Take the minimum of those four total concentrations.  That is your app_solubility at that pH.
	2.	Use scipy.optimize.curve_fit or another optimizer to find the best-fit parameters by minimizing the error between app_solubility(pH_i, params) and the experimental .

4. Example Python Code

Below is a minimal illustrative code snippet using scipy.optimize.curve_fit.  In a real scenario, you will need to:
	•	Decide how you handle  (constant ionic strength? explicit charge balance?).
	•	Possibly refine initial guesses for the parameters.
	•	Add constraints so that all  and  remain positive, etc.

import numpy as np
from scipy.optimize import curve_fit

def total_solubility_neutral(H, S0, pKa):
    """
    If neutral H4A(s) is controlling:
    - [H4A] = S0
    - Then [H3A-], [H2A2-], [HA3-], [A4-] come from acid-base equilibria
    Returns total [A] in solution.
    """
    # pKa = [pKa1, pKa2, pKa3, pKa4]
    pKa1, pKa2, pKa3, pKa4 = pKa
    # Stepwise dissociation constants in "concentration form"
    Ka1 = 10**(-pKa1)
    Ka2 = 10**(-pKa2)
    Ka3 = 10**(-pKa3)
    Ka4 = 10**(-pKa4)

    # We know [H4A] = S0 at saturation
    H4A = S0

    # Then H3A- / H4A = Ka1 / [H+]
    H3A_minus = H4A * (Ka1 / H)
    # H2A2- / H3A- = Ka2 / [H+], etc.
    H2A_2minus = H3A_minus * (Ka2 / H)
    HA_3minus   = H2A_2minus * (Ka3 / H)
    A_4minus    = HA_3minus   * (Ka4 / H)

    return H4A + H3A_minus + H2A_2minus + HA_3minus + A_4minus

def total_solubility_salt(H, Ksp, pKa, charge, M_plus):
    """
    Computes total [A] if a salt of form H_(4-charge)A^(charge-) • (M+)^(charge)
    is the controlling solid.

    For example, if charge=1, we have H3A- controlling, so:
      [H3A-]*[M+] = Ksp  =>  [H3A-] = Ksp / [M+].
    Then use acid-base equilibria to find total [A].

    Input:
    - H = [H+]
    - Ksp = K_sp, e.g. Ksp1
    - pKa = list [pKa1, pKa2, pKa3, pKa4]
    - charge = 1,2,3,... corresponding to how many protons removed from H4A
    - M_plus = [M+] in solution (approx. constant or known)
    """
    pKa1, pKa2, pKa3, pKa4 = pKa
    Ka1 = 10**(-pKa1)
    Ka2 = 10**(-pKa2)
    Ka3 = 10**(-pKa3)
    Ka4 = 10**(-pKa4)

    # The “key” species is H_(4-charge)A^(charge-).
    # We find that concentration from the solubility product:
    #    [H_(4-charge)A^(charge-)] * [M+]^charge = Ksp
    # => [H_(4-charge)A^(charge-)] = Ksp / ([M+]^charge)
    key_conc = Ksp / (M_plus**charge)

    # Now backtrack or forward-track to get all other species.
    # Example: if charge=1, key = [H3A-].
    # From [H3A-], we get [H4A], [H2A2-], [HA3-], [A4-].

    if charge == 1:
        H3A_minus = key_conc
        # H4A
        H4A = H3A_minus * (H/Ka1)
        # H2A2-
        H2A_2minus = H3A_minus * (Ka2 / H)
        # HA3-
        HA_3minus = H2A_2minus * (Ka3 / H)
        # A4-
        A_4minus = HA_3minus * (Ka4 / H)

    elif charge == 2:
        H2A_2minus = key_conc
        # H3A-
        H3A_minus = H2A_2minus * (H/Ka2)
        # H4A
        H4A = H3A_minus * (H/Ka1)
        # HA3-
        HA_3minus = H2A_2minus * (Ka3 / H)
        # A4minus
        A_4minus = HA_3minus * (Ka4 / H)

    elif charge == 3:
        HA_3minus = key_conc
        # H2A2-
        H2A_2minus = HA_3minus * (H/Ka3)
        # H3A-
        H3A_minus = H2A_2minus * (H/Ka2)
        # H4A
        H4A = H3A_minus * (H/Ka1)
        # A4-
        A_4minus = HA_3minus * (Ka4 / H)

    else:
        # If charge=0 or 4, etc., adapt similarly
        raise ValueError("charge not in {1,2,3} in this example.")

    return H4A + H3A_minus + H2A_2minus + HA_3minus + A_4minus

def app_solubility(
    pH,
    S0,
    Ksp1, Ksp2, Ksp3,
    pKa1, pKa2, pKa3, pKa4,
    M_plus
):
    """
    Main function to return the apparent solubility at a given pH,
    given parameters.
    """
    # Convert pH -> [H+]
    H = 10**(-pH)
    pKa = [pKa1, pKa2, pKa3, pKa4]

    # 1) Neutral-limited solubility
    S_neutral = total_solubility_neutral(H, S0, pKa)

    # 2) 1:1 salt-limited
    S_1to1 = total_solubility_salt(H, Ksp1, pKa, charge=1, M_plus=M_plus)

    # 3) 1:2 salt-limited
    S_1to2 = total_solubility_salt(H, Ksp2, pKa, charge=2, M_plus=M_plus)

    # 4) 1:3 salt-limited
    S_1to3 = total_solubility_salt(H, Ksp3, pKa, charge=3, M_plus=M_plus)

    # Actual solubility is the minimum across these possible precipitations
    return np.minimum.reduce([S_neutral, S_1to1, S_1to2, S_1to3])

def fit_solubility_curve(pH_data, S_data, initial_guess, M_plus_const):
    """
    Fits the (pH vs. S) data using scipy's curve_fit.
    :param pH_data: array-like of pH values
    :param S_data:  array-like of experimentally measured solubilities
    :param initial_guess: tuple/list of initial guesses for
        (S0, Ksp1, Ksp2, Ksp3, pKa1, pKa2, pKa3, pKa4)
    :param M_plus_const: chosen [M+] for the fit (or you can incorporate it in parameters)
    """

    # Define a "wrapper" so that curve_fit sees an array for x (pH) and returns model S
    def model_func(pH_array, S0, Ksp1, Ksp2, Ksp3, pKa1, pKa2, pKa3, pKa4):
        # Vectorized evaluation
        return np.array([
            app_solubility(
                pH_val,
                S0, Ksp1, Ksp2, Ksp3,
                pKa1, pKa2, pKa3, pKa4,
                M_plus_const
            )
            for pH_val in pH_array
        ])

    popt, pcov = curve_fit(
        f=model_func,
        xdata=pH_data,
        ydata=S_data,
        p0=initial_guess
    )
    return popt, pcov

# ---------------------------------------------------------------------------
# Example usage:

if __name__ == "__main__":

    # Suppose we have some mock data:
    pH_exp = np.array([2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0])
    S_exp  = np.array([0.1, 0.12, 0.18, 0.25, 0.4, 1.2, 2.5])  # hypothetical

    # We'll guess [M+] is near 0.1 M, for example:
    M_plus_guess = 0.1

    # Initial guesses for parameters:
    # (S0, Ksp1, Ksp2, Ksp3, pKa1, pKa2, pKa3, pKa4)
    init_guess = (0.1, 1e-5, 1e-8, 1e-10, 3.0, 5.0, 7.0, 9.0)

    # Fit
    popt, pcov = fit_solubility_curve(pH_exp, S_exp, init_guess, M_plus_guess)

    print("Fitted parameters:")
    print("S0      =", popt[0])
    print("Ksp1    =", popt[1])
    print("Ksp2    =", popt[2])
    print("Ksp3    =", popt[3])
    print("pKa1    =", popt[4])
    print("pKa2    =", popt[5])
    print("pKa3    =", popt[6])
    print("pKa4    =", popt[7])

    # Then you can plot or otherwise analyze:
    import matplotlib.pyplot as plt

    pH_grid = np.linspace(2, 8, 100)
    S_fit = [app_solubility(
                pH_val,
                *popt,
                M_plus_guess
             ) for pH_val in pH_grid]

    plt.plot(pH_grid, S_fit, label='Fitted curve')
    plt.scatter(pH_exp, S_exp, color='red', label='Experimental data')
    plt.yscale('log')
    plt.xlabel('pH')
    plt.ylabel('Solubility [mol/L] (log scale)')
    plt.legend()
    plt.show()

Notes on the Code and Approach
	1.	Piecewise Minima: The app_solubility function takes the minimum of the solubility limited by each potentially precipitating phase. In practice, one often sees “segment switching” in a solubility–pH plot.
	2.	 Dependence:
	•	In real titration data,  may vary with pH.  If you have a well-buffered system or a known dosing schedule, you can explicitly compute  from the mass balance of added base minus any precipitation.  Or you might treat  as nearly constant at each pH step.
	3.	Bounds & Constraints: In real fitting, you should place constraints so pKa’s remain in ascending order and so that  and  are > 0.  In scipy.optimize.curve_fit, you can specify bounds=([ ... ], [ ... ]) for lower and upper bounds.

5. Summary

By:
	•	Defining each possible solid-phase equilibrium condition,
	•	Converting that to total dissolved species concentrations via the acid dissociation constants, and
	•	Taking the minimum across all possible solids,

one obtains the apparent solubility vs. pH profile.  Fitting that theoretical curve to experimental data (via non-linear least squares) allows one to refine or deduce the unknown parameters .

This general method extends to any polyprotic compound with multiple salts, as long as you include each relevant precipitating phase in the “competition.”